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ABSTRACT 

We present an implicit finite difference representation for general relativistic radiation hydrodynamics 
in spherical symmetry. Our code, AGILE-boltztran, solves the Boltzmann transport equation for the 
■ angular and spectral neutrino distribution functions in self-consistent simulations of stellar core col- 

' lapse and postbounce evolution. It implements a dynamically adaptive grid in comoving coordinates. 

, A comoving frame in the momentum phase space facilitates the evaluation and tabulation of ncutrino- 

^SJ ' matter interaction cross sections, but produces a multitude of observer corrections in the transport 

^ I equation. Most macroscopically interesting physical quantities are defined by expectation values of 

O , the distribution function. We optimize the finite diff'erencing of the microscopic transport equation 

for a consistent evolution of important expectation values. We test our code in simulations launched 
from progenitor stars with 13 solar masses and 40 solar masses. Half a second after core collapse and 
bounce, the protoneutron star in the latter case reaches its maximum mass and collapses further to 
form a black hole. When the hydrostatic gravitational contraction sets in, we find a transient increase 
in electron flavor neutrino luminosities due to a change in the accretion rate. The fi- and r-neutrino 
. luminosities and rms energies, however, continue to rise because previously shock-heated material with 

\^ • a non-degenerate electron gas starts to replace the cool degenerate material at their production site. 

I We demonstrate this by supplementing the concept of neutrinospheres with a more detailed statistical 

, description of the origin of escaping neutrinos. Adhering to our tradition, we compare the evolution 

' of the 13 solar mass progenitor star to corresponding simulations with the multi-group flux-limited 

I diffusion approximation, based on a recently developed flux limiter. We find similar results in the 

. postbounce phase and validate this mgfld approach for the spherically symmetric case with standard 

' input physics. 

, Subject headings: supernovae: general — neutrinos — radiative transfer — hydrodynamics — relativity — 

^ i' methods: numerical 

o : 

^ . 1. INTRODUCTION 

, A s upernova explosion is a dramatic event which includes such a rich diversity of physics l|Bethe I 1199(1 
1^ ' iBurrow s & Young 2000: Mezzacapoa & Brucnn 2000; Janka. Kifonidis & Rampp 2001) that self-consistent numer- 
. !^ ' ical simulations on current computer hardware do not allow to include all relevant pieces at once. After stellar core 
, collapse, a nascent neutron star is formed at the center of the event, requiring a description in general relativity. 
^ ■ Neutrinos generated in the compactifying region are tightly coupled to the matter at high densities, while they leak 
Ci ' or stream out at lower densities. The collapsing stellar core bounces as the equation of state stiffens when nuclear 
densities are exceeded, and a shock wave is formed that ploughs outwards through the still infalling outer layers. Multi- 
frequency radiation hydrodynamics must be used to quantify the energy that the neutrinos deposit in the material 
behind th e shock. This energy deposition has been considered to be essen tial for the success or failure of the supernova 
explosion ijColgate fc White1ll966tlWilson lIToSStlBethe fc Wilson Ill985j) . Observational and theoretical evidence sug- 
gests that this neutrino heating drives convection behind the shock. Instabilities in the protoneutron star, significant 
rotation, and strong magnetic fields further complicate the picture. Observations of neutron star kicks, mixing of 
nuclear specie s, inhomog encous ejec ta, and p olarization of spectra support the presence of as ymmetries in supernova 
explosions (,Tueller et al. 1991: Strom e t al. lfl995: Galama et al~ll 19981 (Leonard et al. II2000I). Motivated by such ob- 
servat ions, the neutrino driven e xplosion mechanism has be en explored in multidimensions Mejant. Benz. fc Co lgatel 
11^92^ 'Mil ler. Wilson, fc Mavle1fl993i iHerant et al. 1 119941: IBurrows. Haves, fc Frvxefl I Il995t l.lanka fc Miiller "19961: 



12002). and alte rnative jet-bas ed explosion scenarios 
Khokhlov et al.1[T999tlMacFadven fc Wooslev Ill99'9l: 



on. fe Mavie I llH9;ii iMerant et al. I lljurro- 
1 119981 iFrver fc Heger II2000I: iFrver fc Warren ll' 



iMezzacap pa et al IFrver fc Heger 

have received n ew momentum l|H6flich. Wheeler, fc Wang I[r999l: 
IWheeler et al. ..2000). 

Because of their excessive computational demand, multi-dimensional simulations have had to rely on physically 
significant approximations. Simulations are performed that impose the neutrino radiation field externally or uni- 
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formly, that prescribe the neutrino spectrum ab initio, that artificially seed instabilities, or that insert jets in order 
to explore important phenomena of supernova explosions in observational data. A complementary approach is to 
model super novae in spherical symmetry. Spherically symmetric simulations of stella r core collapse reach back to the 
late sixties IColga te & White I 1966: Mav & White 1967; Arnctt 1967; Sc hwartz i ri967: Wilson i ITqtH) when com- 
puters first became available. Although the assumption of pure spherical symmetry is also a physically significant 
simplification, it can provide an essential contribution to the quality of numerical supernova models and their inter- 
pretation. A spherically symmetric model can implement full general relativity. Additionally, the reduced number of 
computational zones allows an accurate treatment of neutrino radiation transport and the inclusion of sophisticated 
microscopic input physics. During the search for a robust supernova, mech anism, standards in the nuclear input physics 
l|Lattimer fc Swestv 111991]) and weak interaction physics l|Bruen n 1985) have been established. In recognition of the 
importance of neutrino transport in the supernova, its numerical treatment h as been improved from simple leak- 
age schemes l|Van Riper fc Lat timer lll98lHBaron. Coooerstein. fc Kahana 1ll985t) to multi-group fl ux- limited diffusion 
(mgfld) approximations l|Arnett 11197'^ iB owers fc Wilson 111982^ 'Bruenn lll 985HMvra et ailll987D to investig ations of 
the full Boltzmann transport equation ([Me zzacapua & Brucnn 1993a: Mcss er et al. lll998HYamada. Janka. fc Suzukf] 
urrows et al. 1 12000': Ra mpp fc Janka 112002 ') . Only recently have self-contained simulations of stellar core col- 
lapse and pos tbounce evolution with Boltzmann neutrino tran s port been performed, in Newtonian gravity and the 
0{vlc) limit llRampp fc Janka 11200(1 iMezzacappa et al. ll200lt iThompson. Burrows, fc Pinto Il2003j) and in general 
relativistic space-time l)Liebend6rferet^LTl20 oT'r They allow, on the one hand, the exploration of nuclear physics 
under extreme conditions, an d, on the othe r hand, provide a solid point of reference for the construction of accurate 
multidimensional simulations l|Buras et al. ll2003bt) . 

Basic numerical techniques for radiative t ransfe r have been designed a nd continuously improved over the last three 
decades (see e.g. (jMihalas fc Mihalas1ll984|) and ijLewis fc Miller Ill984|) for a review). In most apphcations, photons 
or neutrons are the transported particles. In our application to supernova dynamics, the transport of neutrinos requires 
a new combination of capabilities derived from both fields, photon and neutron transport. 

As in neutron transport the radiation particles are fermions. In contrast to photons, the neutrinos assume a Fermi- 
Dirac distribution in thermal equilibrium. Weak interactions with nuclear matter can be suppressed due to Pauli 
blocking of final states. This makes the dependence of the collision integral on the distribution function nonlinear. 
Several quantities have to be evolved with special care: An important constraint for the evolution of neutrino and 
electron abundances is lepton number conservation. While applications to photon transport do not require an exact 
count of photons as particles, the outcome of a supernova simulation sensitively depends on the deleptonization. The 
transfer of lepton number occurs via the reactions e~ + ^ n + i^e and e'^ + n ^ p + D^. where the protons, p, and 
neutrons, n, are free or bound in nuclei. The electron neutrino, Ve, and antineutrino, Ve, may escape to infinity or be 
absorbed at distant locations. These reactions determine the electron fraction, and hence, the partial pressure of the 
electron gas in the fluid. This pressure contributes importantly to the dynamics if the electron gas is degenerate and 
the density lower than nuclear density. Beside of the transport of lepton number, the transport of energy is a crucial 
phenomenon in supernova dynamics. Neutrinos escaping from the accreting material and the protoneutron star may 
be absorbed behind the stalled shock. This energy transfer is believed to add to the thermal p ressure behind the 
shock such that ultimately a supernova explosion is launched with the ejection of the outer layers ijColeate fc White! 
ll966tlBethe fc Wilson 111985]) . Last, but not least, the evolution of the global energy should be monitored. Just before 
collapse, the progenitor star is marginally bound. During collapse, the binding energy of the nascent protoneutron 
star increases dramatically. It is balanced by the internal energy of the compressed matter and the trapped neutrinos. 
Only on a longer time scale is the energy of the radiation fleld transported away from the star while the emptied 
states in the neutrino phase space are immediately replenished by neutrino emission in the cooling protoneutron star. 
Since the gravitational energy, the internal energy, and the energy in the radiation fleld rise to the order of 10^^ erg, 
a detailed analysis of global energy conservation is advised in order to make the comparatively small explosion energy 
of lO^^erg predictable. Whether we choose a kinetic equation for the propagation of particle number or an equation 
prescribing the evolution of radiation intensity, it is a numerical challenge to conserve both lepton number and energy 
with the same finite difference representation of the transport equation. 

The following transport requirements more closely resemble the astrophysical applications of photon transport than 
neutron transport. Owing to the protoneutron star or black hole formed at the center of the event, space-time is 
curved. The particles follow geodesies in general relativistic space-time, i.e. their angular distribution is affected by 
gravitational bending. Additionally, the particle energies are subject to a gravitational frequency shift. In contrast 
to many neutron transport applications, the fluid is highly dynamic during collapse and even more so after shock 
formation. Thus, observer corrections for Doppler shift and angular aberration have to be implemented as in photon 
transport applications. However, there is some freedom in the choice of where to apply these corrections. In an inertial 
frame they complicate the description of interactions in the collision integral. In a comoving frame they enter the 
transport terms in the Boltzmann equation. 

We proceed with the latter choice in an S'Ar-method. Our code, AGILE-boltztran, emerged from the following com- 
ponents: The neutrino transport part, BOLTZTRA N, has been developed for the sim ulation of stellar core collapse in an 
implicitly flnite differenced 0{v/c) approximation (jMezzacappa fc Bruenn ll993albO . It has been compared to mgfld 
in selected stationary state phases, and standard test problems for radiative transfer in supernovae have been performed 
l|Messer et al. f 1998: Messcr 2000). AGILE is an implicit general relativistic hydrodynamics code t hat evolve s the Ein- 
stein equations based on conservative finite differencing on an adaptive grid (Liebendorfer. Rosswog fc Thic lcm anrTI 
I2002D . In this paper, we describe how these codes are merged and extended to enable accurate simulations of the very 
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dynamic postbounce phase. In particular, we detail the finite differencing of the observer corrections for the dynamic 
conservation of particle num ber and total energy in the transport scheme, and describe the extension of BOLTZTRAN 
to general relativistic flows l)Liebend5rfer ||200C1|) . In section |21 we start with the characterization of the equation of 
state and a list of the included neutrino-matter interactions, followed by a collection of the basic equations of general 
relativistic radiation hydrodynamics in spherical symmetry. Two exemplary runs from collapse through bounce and 
postbounce evolution are also described in this section. In or der to test our code in a bro ad range of conditions, 
we launched simulat ions from a small 13 MfD progenitor star l)Nomoto fc Hashimoto l[l988() and a very massive 40 
Mq progenitor star ijWooslev fc Weaver Ill995|) . We investigate the different regions of neutrino emission in the star. 
Section O is entirely devoted to the documentation of the finite differencing in our computer model. In section 01 
we analyze and verify the performance of our code in various example situations encountered in the evolution of the 
two simulation runs. We close the section with the comparison of our results with Br uenn's multi-group flux-limited 
diffusion code which implements a recently developed new flux limiter l|Bruenn 112001 . 

2. PHYSICS IN THE MODEL 

The physical model can be divided into two parts. On the one hand, there is the microscopic physics input-the 
specification of particle abundances and reaction cross sections. The microphysics in supernova models is continuously 
improving and many uncertainties remain to be resolved. However, in this methodological paper, we will only shortly 
summarize the ingredients that were standard at the time the code was written and hope that our code AGILE- 
BOLTZTRAN will Continue to be useful in future discussions and evaluations of input physics improvements. On the 
other hand, there are the radiation hydrodynamics equations we are solving on the computer. They are the foundation 
of our implementation of neutrino transport and receive a detailed discussion later in this section. 

2.1. Equation of state and weak interactions 

The equation of state describes the thermodynamical state of a fluid element ba sed on density, p, temperature, 
r, and the composition. We use the equation of state of lLattimer fc Swestvl l)1991|) . It assumes nuclear statistical 
equilibrium and we apply it wherever the density is larger than 10^ g/cm"^ and the temperature larger than 5 x lO^K. 
This region is described by a liquid drop model for a representative nucleus with atomic number A and charge Z, 
surrounded by free alpha particles, protons, and neutrons. The baryons are immersed in an electron and positron 
gas in equilibrium with a photon gas. Beyond nuclear density , where no isolated nuclei are present, the complicated 
population of hadrons l|01endenning IIT985I : iPons et al. Ill999j) is approximated by bulk nuclear matter comprised of 
protons, neutrons, and electrons. However, the central density of the protoneutron star at bounce reaches only about 
twice nuclear density and the hadron population may only develop later, after the very dynamical postbounce phases. 
At the low temperature border of nuclear statistical equilibrium, the equation of state is connected to a Boltzmann gas 
of silicon atoms. In any of these cases, once the density and temperature are given, the composition is fully determined 
by the specification of the electron fraction Yg- 

Matter is connected to the neutrino radiation field by weak interactions. We consider neutrinos of all three fiavors 
and assume that they are massless. The weak interactions enter the collision term in the Boltzmann equ ation as energy - 
and angle-dependent emissivities, opacities, and scattering kernels. We include the set specified by l)Bruenn II198,'t() : 
(i) electron-type neutrino absorption on neutrons, (ii) electron antineutrino absorption on protons, (iii) electron-type 
neutrino absorption on nuclei, (iv) neutrino-nucleon scattering, (v) coherent scattering of neutrinos on nuclei, (vi) 
neutrino-electron scattering, and (vii) neutrino production from e lectron /positron pair an nihilation. These re actions, 
and their inverses, are implemented in our code as described bv iMezzacappa fc Bruenn" ( 1993b d): iMesserl |2000). 
In the following code description, we will only include emissivities, j, and opacities, x? because this is sufficient to 
describe how the collision term enters the transport and hydrodynamics equations. In the simulations, the scattering 
kernels are included in the collision integral as well. 

The particles treated by the equation of state are assumed to react with each other on very short time scales such 
that a description in terms of an instantaneous equilibrium is appropriate. Neutrinos in high-density regimes can also 
achieve local thermal and weak equilibrium with matter if the opacities are sufficiently high. Unlike the equilibrium 
with respect to the strong interaction, however, this equilibrium must be determined within our solution of the 
transport equation. For example, in the protoneutron star at densities above 10"'^^ g cm~'^ and temperatures above 
5 X lO^*' K the neutrinos are trapped and are well-described by a Fermi-gas in thermal equilibrium with the fiuid. At 
lower densities, the thermalization time scale becomes longer; then the neutrinos can propagate with a nonequilibrium 
spectrum throughout these regions, to be absorbed elsewhere or leave the star. The strong coupling of the neutrinos 
to the matter at high densities and the strong coupling between different locations mediated by neutrino transport 
complicates the evolution of a numerical solution. If the problem is separated into independently updated pieces by 
operator splitting, the numerical solution will only be stable if information in the numerical implementation is shared 
faster between the independent updates than in the evaluated physical processes. The fast time scale of neutrino-matter 
interactions and the propagation of neutrinos at light speed may severely restrict the time step. The required coupling 
can be built directly into the numerical scheme by an implicit finite differencing of essential parts of the transport 
equation. Unfortunately, such a differencing requires knowledge of the derivatives of the collision term with respect to all 
independent state variables-i.e., density, temperature, electron fraction, and neutrino distribution functions. Because 
the emissivities, opacities, and scattering kernels strongly depend on the neutrino energies and, in the scattering case, 
on the neutrino propagation directions, the numerical evaluation of the collision term and its derivatives becomes a 
nonnegligible part of the overall computational effort. IMezzacappa fc Bruenn I l)1993a(l developed a storage scheme 
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that allows the reuse of previously calculated emissivities, opacities, and scattering kernels by linear interpolation 
within a dynamical table in the independent variables of logarithmic density, logj^g logarithmic temperature, 
logj^Q (r), and electron fraction, Y^. If one uses these same independent variables in the implicit formulation of the 
Boltzmann equation, the correct partial derivatives of the reactions directly emerge from the coefficients of the linear 
interpolation, without additional computational effort. On the one hand, the reuse of previously evaluated interactions 
is straightforward if the transport equation is solved in the rest frame of the fiuid, such that no transformation of the 
neutrino energy or angle dependence of the interactions is required. On the other hand, the transport equations are 
simpler in the laboratory frame. In this paper we demonstrate that in the case of spherical symmetry the complexity 
of the transport equation is manageable and proceed with the analysis in a comoving frame (spacetime coordinates 
and neutrino four-momentum) to take advantage of the simplifications in the collision term. On average, we have to 
evaluate new collision integrals in about two to three zones per time step (out of a hundred zones). Although these 
numbers depend very much on the specific phase of the simulation, the evaluation of nuclear physics input may still 
take about half of the total execution time. 

2.2. Radiation hydrodynamics in spherical symmetry 

Many sphericaUy symnietric simulati ons of compact obj ects have been approached in comoving orthogo- 
nal coordinates (iMisner & Sharo 1964 jM^y ^hite | [19 6^. Finite d i fference schemes of vary i ng complexit; 



were designed in (Mav & White 1967; Van Ri per 1 119791: IBruenn 1 119851: iR.ezzoll a fc Miller 1 119 94: 'Swestv E 

llYai 



ILie bendorfer. Rosswog & Thielemann 2002), culminating in an approximate Riemann solver ijYamada ..1997^ . The 
left-hand side of the Einstein field equation, the Einstein tensor, is based on the metric 

ds^ = -a^dt^ + (^Y^ da^ + {dd"^ + sin^ dd^p^) , (1) 

where r is the areal radius and a is a label corresponding to an enclosed rest mass (the prime denotes a derivative 
with respect to a: r' = dr/da). The proper time lapse of a comoving observer is related to the coordinate time dt 
by the lapse function a. We have made the substitution gaa = r' /T, based on a function F {t, a), for the space-space 
component of the metric. The angles i) and (p describe a two-sphere. We use natural units such that the velocity of 
light, c, and the gravitational constant, G, become 1. 

The right-hand side of the Einstein equations is given by the fluid- and radiation stress-energy tensor, T. In a 
comoving orthonormal basis, it has the components (Lindquist 1966) 

T" = p{l + e + J) 
T'"'=p + pK 

T^^ =T'P'^=p+^piJ - K). (2) 

The total energy is expressed in terms of the rest mass density, p, the specific internal fiuid energy, e, and the specific 
radiation energy, J. The isotropic fiuid pressure is denoted by p, and the radiation stress is composed from the zeroth 
(J) and second (K) angular moments of the radiation intensity. Radial net energy transport is accounted for by the 
nondiagonal component of the stress-energy tensor, the first angular moment (H) of the specific radiation intensity. 

We define a velocity m, equivalent to the r component of the fluid four- velocity as observed from a frame at constant 
areal radius r ijMav fc White iri967tl . and identify the total energy enclosed in a sphere with the gravitational mass, 
m. In the special relativistic limit, F = y^l -\- u'^ — 2m /r then becomes the Lorentz factor corresponding to the boost 
between inertial and comoving observers. As in nonrelativistic hydrodynamics we can define a specific volume, 1/-D, 
specific energy, r, and specific radial momentum, 5, by 



(3) 

1 



1 _F 

r^Tie + J) + ^i^-u'-fj+uH (4) 

S = u{l + e + J) + rH. (5) 

It has been shown in ijLiebendorfer. Mezzacappa. fc Thielemannl 12001!) that these definitions lead to conservation 
equations iQ-® that are analogous to the continuity equation, the conservation of total energy, and the conservation 
of radial momentum^ : 

This reference unnecessarily assumes an isotropic radiation stress. However, we note that the difference between the full stress-energy 
tensor and the isotropic approximation has exactly the same form as the artificial viscosity tensor introduced in the same reference to 
numerically stabilize shock fronts. Hence, in all derivations in the above reference we may simply use the pressure p = p + pJ/3 for the 
isotropic part and set the viscosity coefficient, Q, to Q = — p(J/3 — K), in order to obtain a description of radiation hydrodynamics that 
extends to the case where large radiation energies do not satisfy J 3K. 
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The change of the specific volume in Eq. @ is given by the balance in the displacement of the zone boundaries. The 
rate of change of total energy in Eq. Q is determined by the surface luminosity, L = Anr'^ pH , and the work on 
the surface of the mass shell against the pressure, p + pK. Of leading order in the momentum equation (jSJ are the 
pressure gradient and the gravitational force, m/r^ . The constraints Q and are most easily understood in the 
Newtonian limit (the enclosed volume is defined hy V = 47rr'^/3), where the first becomes the definition of the rest 
mass density and the second the Poisson equation for the gravitational potential. The time derivative in equation Hll|l 
is very small; therefore, this equation essentially acts as a constraint on the lapse function, a. This equation derives 
from the space component of the four-divergence of the stress-energy tensor. In addition to the evolution of the total 
energy, we also need an equation for the evolution of the internal energy that we may derive from the time component 
of the four-divergence of the stress-energy tensor: 



d_ 

dt 



[e J] = [ATTT^a^pH] - {p + pK) 



d_ 

dt 



au 



(J - 3if) , 



(12) 



Next, we detail the description of the radiation field. We identify the energy flux pH with a particle flux that is 
determined by a Boltzmann transport equation. The transport equation is split into a left-hand side and a right-hand 
side. The left-hand side is the directional derivative of the particle distribution function along trajectories of free particle 
propagation. This derivative is equated to the changes in the distribution function du e to co llisions, which are described 
by the right hand sid e of the equation. Once a 1-1-1 decomposition of space-time ijArnowitt. Deser fc Misner1ll962t 
iSmarr fc Yorkiri978|) and a basis in the momentum phase space for the particle four-momentum have been chosen, 
the directional derivative along the phase flow can be expressed in terms of partial derivatives of the distribution 
function with respect to the space-time coordinates and momenta (Lindquist 1966; Mczzac appa fc Matzncr 19891. 
We measure the particle four-momentum in a comoving orthonormal frame, with components 



pcosu, p 



i5 



p sin 9 COS (I), p'^ = p sin 9 sin (j). (13) 

In spherical symmetry, the particle energy, E, measured in a comoving frame, and the cosine of the angle between the 
particle momentum and the radial direction, p = cos 9, completely describe the particle phase space. The neutrinos are 
assumed to have no mass. In spherical symmetry, the distribution function does not depend on the three-momentum 
azimuth angle (p. Thus, the specific particle distribution function depends on four arguments and describes the number 
of particles at a given time, t, in the phase space volume E'^dEdpda by 

F{t,a,p,E)E'^dEdpda. (14) 

Boltzmann equation reads l)Yamada. Janka. fc SuzukTI Il999t 



With the metric of Eq. 



dN 
the 



ILiebendorfer. Mezzacappa. fc Thielenianrni200H) 



Ct+Da+D.+DE + Ou + OE= Cc, 



(15) 



with 



Ct 



dF 
adt 



a da ■' 



(16) 
(17) 
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Cc = --xF. (22) 
P 

The source on the right-hand side, Cc, is the colhsion term that describes changes in the particle distribution function 
due to local interactions with matter. It is represented here by an emissivity j and an opacity x- AH other terms stem 
from the partial derivatives of the distribution function with respect to the phase-space coordinates in the directional 
derivative along the phase flow. They can all be physically interpreted. The first term on the left hand side of the 
equation, C(, is the temporal change of the particle distribution function. The second term, Da, counts the particles 
that are propagating into or out of an infinitesimal mass shell. The third term, D^, accounts for the change in the 
neutrino distribution function in an angle interval owing to the propagation of the neutrinos along geodesies with 
changing local angle cosine /x. The curved particle trajectories in general relativity are accounted for by the term 
proportional to the gradient of the gravitational potential, $, 

Ida _d<^ 
a dr dr 

The fourth term, De, expresses the redshift or blueshift of the particle energy that applies when the particles have a 
velocity component in the radial direction (/i ^ 0) and, therefore, change their position in the gravitational well. The 
fifth and sixth term, Oe and O^t, account for the Doppler shift and the angular aberration between adjacent comoving 
observers. 

The integration of the Boltzmann equation over momentum space, spanned by the particle direction cosine and 
energy, gives the local conservation laws for particle number and energy. We define and to represent the 
zeroth and first /j, moments of the distribution function: 

1*1 poo 



1 Jo 

pi />oo 



H =J J FE'^dEfidfj.. (23) 



^ Uirr^apH^] ~a [ ^E^dEd^ + a [ xFE^dEdfi = 0. (24) 
da ^ ' J p J 



Integration of Eq. (|15|l over p and E with E^ as the measure of integration gives the following evolution equation for 

dt 

The derivatives with respect to the momentum phase space in Eq. \\b\ do not contribute because (l — /i^) vanishes 
at /Lt = ±1 and E^F is zero for = and E — oo. Eq. H24|l is a continuity equation analogous to Eq. ((HJ, extended 
by source and sink terms for the radiation particles. One more integration over the rest mass a from the center of the 
star to its surface gives the evolution equation of the total particle number. 

Slightly less straightforward is the derivation of total radiation energy conservation. We define the energy moments 





f FE^dEdp, 


"-J 


[ FE^dEfidfi 




[ FE^dEp^djjL 




[ FE^dEfi^dfi 



(25) 

and evaluate the evolution of the radiation energy as measured by an observer at infinity, 

^ J (T + ufi)FE^dEdn. (26) 

To this purpose, we integrate Eq. (|15|) again over phase space, but this time with measure of integration (F + up.) E^. 
After performing some integrations by parts to account for the time and space dependence of F and u, this leads to 
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the concise result ijLiebendorfer. Mezzacappa. fc Thielemann 11200 Ij) 



d ( ^ 

= — (r J + uH) + — {A.Txr'^ap {uK + VB)\ + 47rrap ( 1 + e + - ) 

-ar / f ^ - X ) E^dEdn + au I xFE'^dEfidfi. (27) 



Note that the conserved quantity TJ + uH is the radiation energy density in the frame of an observer at infinity. It 
is expressed in terms of the momentum moments J and H in the comoving frame. The second term describes the 
surface work by the radiation pressure, pK , and the energy loss or gain due to the luminosity L = Anr'^pH at the 
boundary. The third term contains a gravitational term coupling the matter enthalpy with the luminosity that we 
neglected in previous work (Lie bendorfer. Mezzacappa. fc Thielemann 11200 111 . The source terms in Eq. H27|l describe 
the energy exchange with matter by particle emission, absorption, and radiation stress. The omitted terms from 
neutrino scattering enter the equation in a similar form. Beforehand, we found that Eq. ITJ describes the evolution of 
the total energy. Then, from the Boltzmann equation, we derived Eq. H27(l for the evolution of the radiation energy. 
Thus, we will find a consistent equation for the evolution of the hydrodynamics part by the subtraction of Eq. H27|l 
from Eq. Q. The result is Eq. H89|) in section [3 . 41 where we discuss the implementation of the hydrodynamics part. 
The same procedure applied to other conserved quantities leads the full set of consistent hydrodynamics equations. 
By taking moments of the Boltzmann equation with the measures of integration E^ {u + Tp), E^p/ (^Anr^p), and E^ , 
respectively, we derive equations for the evolution of the radiation momentum, Eq. 1)28(1 . a radiative contribution to 
the lapse function, Eq. H29|l . and the radiation energy in the comoving frame, Eq. H3(J|I : 

— [uJ + rH]^-— [iirr^ap {TK + uH)] 

"r(j + 3if)™-fl-^') iJ-3K)-2K 
r \ r J 

+ 47rr2 (J {p + pK) + {\ + e + J) pK - 2pH^)] 

-aV j xFE^dEpdp + au J (J- - xF^ E^dEdp (28) 



d 



dt l^-Kr^p 



H 



oa p oa iVD 



1 



xFE'^dEpdp (29) 



- = ^-[,.rWpH]-pK- [-y-iJ^m 

+ aj - xF^ E^dEdp. (30) 

The subtraction of Eq. (|28|) from Eq. ||SJ| leads to a hydrodynamics equation for the evolution of the momentum, 
Eq. H90|l . The subtraction of Eq. (|29|) from Eq. Hll|) leads to a hydrodynamics equation for the update of the lapse 
function, Eq. (|94() . Finally, the substraction of Eq. 130() from Eq. (|12ll leads to a hydrodynamics equation for the 
evolution of the internal energy, Eq. ((95|l . We will pay attention to preserve this consistency also in our finite difference 
representation of the equations of radiation hydrodynamics. However, before we proceed with the technical details in 
section 121 we provide in the next subsection an overview of two exemplary simulation runs to complete the physical 
context and to illustrate the numerical challenges we face. 

2.3. Neutrino transport in two representative simulations 

In this section, we provide an overview of the core collapse and postbounce evolution in our models for the 13 
Mq and 40 Mq stellar progenitors. These provide the physical context for t he code tests in section IH A tho rough 
discussion of supernova physics can be found in previou s reviews, e.g. in llBethelll99(Tt H urrows fc Young] 12000; 
iMezzaca ppa fc Bruenn ll200CtlJanka. Kifonidis fc Rampp l HoOl). Earlier runs of the 13 Mq model have been described 
in f Mezzacapp a et al. 11200 ll : iLiebendorfer et al. I with Newtonian and general relativistic gravity. Here, we add 

information on the formation of the neutrino spectra and report on our first self-consistent simulation running all 
the way from core collapse to the onset of black hole formation in the case of the 40 M© model. We have chosen 
progenitors on the light and massive side with respect to the range of potential core collapse supernova progenitors. 
This demonstr ates the spread in the results. The results of o ur simulations for intermediate mass progenitors are 
summarized in l|Messer Il200"ol ILiebendorfer et al. ll2001all2002D . The latest was calculated with the finite differencing 
described in this paper. 

The most prominent characterization of a supernova explosion is the trajectory of the shock position. We define 
the shock position as the location with the maximum infall velocity (i.e. minimum in the velocity profile). Before the 
shock is formed after bounce, the location with maximum infall velocity coincides with the sonic point which separates 
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Fig. 1. — Shown are the shock trajectories of the 13 Mq (dashed hne) and 40 Mq (sohd line) models. At negative times, i.e. before 
bounce, the trajectories indicate the position of the sonic point instead of the position of the not yet formed shock wave. The left hand 
part of the figure zooms in on the time around bounce to demonstrate that the shock is formed at the sonic point (no discontinuity in the 
lines) and that this happens at the same position in both models. The right hand part shows the shock trajectories over a longer time 
scale. After about 500 ms, the 40 Mq star collapses to a black hole. The perturbation at tpj, =~ 0.3 s is the consequence of a physical 
change in the luminosities. The two perturbations after ipj, =~ 0.4 s are the result of an artificial increase of the numerical shock width 
we had to apply in order to run the simulations to the end. 



the causally connected inner core from the supersonically infalling outer core. The pressure wave emerging from the 
center at bounce turns at this transition point into a shock wave. The trajectory of maximum infall velocity is therefore 
continuous across bounce as shown in the magnifying window on the left hand side of Fig. If we compare the 

position of the sonic point in the 13 Mq model with its position in the 40 Mq model, we find that they converge to the 
same point before bounce. For the explanation, we recall that the sonic point depends on the Chandrasekhar mass, 
which is determined by the electron fraction profile. The electron fraction profil e of the two runs converges due to a 
strong feedback of the electron fraction on the free proton abundance (Mcsscr e t al. 11200,31) . Within the "s tandard" 
input physics, it is assumed that electron capture on nuclei with a full neutron f7/2 shell is Pauli- blocked (jBruenn I 
Il985j) . Under such conditions, our simulations allow only electron capture on free protons. Now, if the electron 
fraction and/or temperature in one model would only be slightly higher than in the other model, this would result in 
a significantly higher free proton abundance. It would cause a significantly higher number of electron captures than in 
the other model. Hence, the differences between the models are reduced. It has recently been shown that, due to the 
finite temperature in the nuclei a nd due to co rrelations, electron capture on nuclei dominates electron capture on free 
protons throughout core collapse l|Langanke et al. 2003). It will be interesting to see if the described feedback will to 
the same extent be at work with these more realistic electron capture rates. We can only expect this if the electron 
capture on low abundance nuclei would turn out to be comparable to, or larger than, the electron capture on the most 
abundant nuclei (the quantity to compare would of course always be the product of the electron capture rate with the 
abundance of the target). The right hand side of Fig. Q shows the shock trajectories over a longer time interval, 
up to one second after bounce. The shocks recede in both models after tpb = 100 ms. The conditions for a shock 
revival deteriorate. The intensive neutrino emission from the cooling region undermines the pressure support below 
the heating region and the material is drained from the latter onto the protoneutron star ijJanka Il2001() . The transient 
stall in the receding shock front in the 40 Mq model at 0.3 s after bounce is a reaction to an enhanced electron fiavor 
luminosity. It will be further analyzed below, together with the description of the evolution of the luminosities. In the 
evolution of the 40 Mq progenitor, we encountered a numerical stability problem after tpf, = 0.4 s. The adaptive grid 
created extremely small mass zones such that the convergence radius of the Newton-Raphson algorithm in the implicit 
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Fig. 2. — The luminosities and rms energies of the neutrinos are shown as a function of time. The results of the 13 Mq model are drawn 
with dashed lines and the results of the 40 Mq model with solid lines. A thick line belongs to the electron neutrino, a line with medium 
width to the electron antineutrino, and a thin line to the fj,- and r-neutrinos. We sampled the luminosities at a radius of 500 km. Also 
the rms energy of the neutrino flux (not abundance) was calculated at this location. We adjusted the time coordinate by At = —500 km/c 
to account for the (approximate) propagation time to the sampling radius. The two progenitors show a comparable neutrino burst with a 
peak height of 3.5 X 10^^ erg/s. Significant differences appear in later phases. The variations in the density profiles in the outer layers of 
the two models determine the accretion-dominated electron flavor luminosities. 



hydrodynamics was severely reduced due to truncation errors. We increased the artificial viscosity in two sequential 
steps to widen the shock and continue the run. This numerical shock widening is responsible for the outward steps in 
the shock position ~ 0.5 s after bounce. Shortly after the collapse of the protoneutron star to a black hole has set in, 
our code crashes unavoidably because of the coordinate singularity in the comoving coordinates at the Schwarzschild 
horizon. We will extensively use the hydrodynamic profiles at 0.4 s after bounce for testing in later subsections. 

The neutrino luminosities and rms energies are shown in Fig. (|2Jl. The electron neutrino luminosities rise during 
core collapse and reach a level of 10^^ erg/s. The collapse is halted at nuclear densities and a bounce shock propagates 
outwards through neutrino opaque material. The neutrino luminosity decays to a 30% lower level during this short 
period of ^ 4 ms duration. It has b een assigned to a decrease in the free proton fraction when the shock is formed 
ijThompson. Burrows, Pin to 20031. Additionally, while the shock is running out to the neutrinospheres, it condenses 
previously still neutrino emitting material to even more neutrino opaque densities. When the shock reaches the electron 
neutrinosphere, an electron neutrino burst with a peak height of 3.5 x 10^'^ erg/s is launched by copious electron capture. 
As the neutrinos escape quickly, the freed phase space is refilled with new neutrinos and the matter deleptonizes rapidly. 
This phase is very similar in both models because, after core collapse, the structure of the inner core is so similar. 
Differences appear later, when the accretion luminosity dominates over the core diffusion luminosity with a ratio of 
about 2 : 1 to 3 : 1. The higher densities in the outer layers of the more massive progenitor produce considerably 
higher accretion luminosities when they settle in the gravitational potential on the surface of the protoneutron star. 
The rising electron neutrino rms energies before bounce reflect the conditions in the compactifying material at infall. 
After the neutrino burst, the rms energies adjust to the spectrum set by the shock heated mantle and reflect the 
conditions at the location of decoupling. However, we have to keep in mind that the location of decoupling strongly 
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varies for individual neutrino flavors and neutrino energies. Thus, the rms energy rather reflects an emission-weighed 
samphng of conditions at different locations. Quite generally, the fi and r neutrinos decouple deeper because of the 
insensitivity of these neutrinos to charged current reactions. The electron antineutrinos decouple deeper than the 
electron neutrinos because of the smaller proton than neutron abundance. 

In order to investigate the origin of the neutrino luminosities in more detail, we introduce in appendix IbI radius- and 
energy-dependent attenuation coefficients, ^ (r, E), that express the probability that a neutrino with energy E emitted 
at radius r escapes from the computational domain without a reaction that changes its type or energy. The attenuation 
coefhcients carry information tha t is similar to the optical depth, T{r,E): £_(r,E) ~ exp (— r (r, i?)). However, their 
evaluation according to Eq. (|B5|) is fully consistent with the finite differencing of the Boltzmann equation and takes 
the variations in the local flux factors into account. Assume for example, we investigate a reaction £ that produces at 
radius r neutrinos in the energy interval dE with an energy emissivity E^j^{r, E)dE. With the help of the attenuation 
coefficients, we can quantify the contribution of a given volume element Airr'^dr to the total luminosity, 



The total neutrino luminosity of a given neutrino type is then given by the integral of (r, E) over energy and position 
for all reactions that contribute. 



We demonstrate in appendix that the attenuation coefficients in Eq. I)B5(I represent the total luminosity in the 
simulation accurately. In the following figures, we cumulatively plot the quantity (r, E) AE in units of erg/s/km. It 
is natural to choose AE in accordance with the width of the fcmax = 12 energy groups we used in the simulations. At 
the bottom of the figure we start with the electron or positron capture reaction (depending on the neutrino type) . We 
draw g^'^P^ (r, Ei) AEi for the lowest energy group and enclose it by a blue line. On top of it, we add g'^^P^ (r, E2) AE2 
for the next energy group, again enclosed by a blue line. The shading of the enclosed areas indicates the energy group 
according to the legend at the bottom of the figure. Energy groups that do not contribute collapse to a single line 
with zero enclosed area. On top of g'^^^*' {r, Ek^_^^) AEk^_^^ we continue with gP'^" {r, Ei) AEi for the pair creation 
reaction. We enclose this and the following area elements by a white line to distinguish them from the electron or 
positron capture reactions. After the addition of gP'^" {r, Ek^^^) AEk^_^^, we continue with the contributions from 
neutrino-electron scattering, i.e. g^^^^ (r, Ei) AEi to gS'^^' (r, E^^^^) AEk^^^- We use green lines as separators for the 
neutrino-electron scattering. The figures become intuitively accessible as soon as one realizes that the total shaded 
area is proportional to the total luminosity of the star. The total area of a specific energy shading is proportional to 
the contribution to the total luminosity of neutrinos from the corresponding energy group. The total area of a specific 
reaction is proportional to the contribution to the total luminosity of neutrinos with this last inelastic interaction 
before the escape. A cross section through the shaded area at a given radius tells about the spectrum of the neutrinos 
escaping from that region, and about the probability of the reaction type they had at that position before the escape. 
In order to characterize also the extent of isoenergetic scattering of neutrinos off nucleons and nuclei, we mark the 
neutrinospheres for the energy groups at the top of the figure. The energies are rising from the left to the right 
according to the legend at the bottom of the figure. For the interpretation of the figures it is also useful to know 
the thermodynamical conditions at the locations the neutrinos are emitted. For each density decade we set a marker 
at the bottom of the figure. Additionally, we include the electron fraction in the graph with the electron neutrino 
analysis, and the entropy in the graph with the electron antineutrino analysis. The solid line represents the profile 
in the simulation, the dashed line represents the equilibrium value that would be achieved by infinitely long exposure 
of the stationary fluid element to the prevailing neutrino abundances. Finally, the graph with the /i- and r-neutrino 
analysis obtains profiles with the temperature (dashed line) and electron chemical potential (solid line). 

Fig. ||2Jl, for example, shows a snapshot at 5 ms before bounce in the collapse of the 13 M0 progenitor. Neutrinos are 
escaping from the range between 20 km and 200 km radius, roughly corresponding to a density range of 10^" g/cm^ to 
10^^ g/cm'^. The electron fraction profile (solid line) reflects the deleptonization that has already occurred in the more 
interior regions where Ye approaches values around 0.3. The equilibrium value (dashed line) is still lower, indicating 
that the deleptonization is ongoing and that the deleptonization time scale is slightly slower than the dynamical time 
scale. The pair process does not contribute in this collapse phase of high electron degeneracy. The corresponding area 
collapses to one white line in the graph that separates the electron capture contributions below it from the neutrino- 
electron scattering contributions above it. Almost no neutrinos escape directly after an electron capture at a radius 
as small as 50 km. Neutrinos escaping from this region are thermalized by scattering off electrons. They escape with 
quite low energies. Around 100 km radius we find that about half of the escaping neutrinos stem directly from electron 
capture, while the other half has scattered off an electron. At 150 km '--^ 2/3 of the neutrinos escape without further 
electron-scattering. Among the neutrinos produced by electron capture, the neutrinos with higher energies escape 
from smaller radii than the neutrinos with lower energies. In this collapse phase, the "standard" input physics only 
includes electron captures on free nucleons. The Q value of the reaction is small and very few low energy neutrinos 
are directly produced if the electron chemical potential is large (e.g. 17 MeV at 100 km radius in this time slice). 
The escaping low energy neutrinos have scattered off electrons. However, the larger Q value of more realistic electron 
capture rates on neutron-rich nuclei may shif t the ene rgy of directly escaping neutrinos to lower values such that the 
count of direct escapes is increased (^Langanke et al. ll2003 ). Finally, we remark that the neutrinos in a given energy 
group are produced at a significantly smaller radius than the location of the corresponding transport sphere. This is 



5^ (r, E) dEdr = ^(r, E)f (r, E) E^dEAirr'^dr. 
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Fig. 3. — The last inelastic interactions of escaping neutrinos at 5 ms before bounce in the 13 Mq model. The abscissa of the graph is 
the radius, ranging from the center of the star to 250 km radius. The density markers at the bottom of the graph indicate the position 
of density decades, logj^Q (p), where p is given in g/cm^. The thick solid line shows the electron fraction profile. The thick dashed line 
shows the electron fraction in weak equilibrium under otherwise unchanged conditions (same density, temperature, neutrino abundances, 
and spectra). These quantities belong to the ordinate at the left hand side. The total shaded area in the graph corresponds to the total 
electron neutrino luminosity at 5 ms before bounce in units of erg/s. The appropriate ordinate on the right hand side then carries the 
units erg/s/km. The shaded area is divided into three sections, according to the type of the last energy-changing reaction of the escaping 
neutrinos. Segments separated by blue lines in the lower part of of the shaded area outline the contribution of electron captures to the 
escaping neutrinos. The contribution from pair production is bordered by white lines. The contribution from neutrino-electron scattering 
is shown in the upper part of the shaded area, bordered by green lines. As there are no contributions from pair production in the cool 
and electron-degenerate material at 5 ms before bounce, the section belonging to the pair production reaction collapses to a white line 
between the electron capture section and the neutrino-electron scattering section. The contribution for each reaction is further subdivided 
into contributions from each energy group in the simulation. The intensity of the shading identifies the neutrino energy according to the 
legend at the bottom of the figure. Each energy group has its own neutrinosphere at optical depth t = 2/3. Their locations, marked in 
the upper part of the figure, statistically indicate the positions of the last interaction of the neutrinos with matter, isoenergetic scattering 
included. The energies rise from the left to the right according to the legend. The transient ondulations in the luminosity contributions 
in this phase are a numerical artef act caused by the low resolution of the Fermi surface in the energy dimension of the momentum phase 
space (see discussion in subsection 14.51 . 



the result of the dominance of the isoenergetic scattering cross section in the collapse phase. The diffusive propagation 
of the neutrinos extends to much lower densities than neutrino-electron scattering or neutrino absorption. 

We switch to the next interesting phase: the electron neutrino burst at '--^ 5 ms after bounce (Fig. While the 
region of neutrino emission during core collapse was very broad, it is extremely narrow (~ 10 km) in the burst phase. 
In the neutrino burst, the electron neutrinos escape directly from electron capture. Some neutrino-electron scattering 
does occur in front of the shock in an earlier st age and inelastic scattering of burst neutrinos on nuclei in front of 
the shock are possible ijBruenn fc Haxtoii1ll99l|) . but not included in the standard input physics. Only the energy 
groups at 10.5 MeV and 16 MeV contribute significantly to the burst. This is in good agreement with the position of 
the corresponding transport sphere in the upper part of the figure (fourth and fifth from the left, respectively). The 
density is of order 10^^ g/cva^ . The trapped electron neutrinos inside the region of main emission are in weak (and 
thermal) equilibrium with the fluid. This is evident in the congruence of the electron fraction profile (solid line) with 
the equilibrium Ye (dashed line). 

The situation at 50 ms after bounce is shown in Fig. This is the phase where neutrino heating starts to set 

in. All luminosities are fully developed at this time. Graph (a) reveals electron capture as the almost exclusive source 
of electron neutrinos. The higher energy neutrinos emerge from larger radii. The isoenergetic scattering cross sections 
are comparable to the neutrino absorption cross sections. Thus, the regions of emission for the different energy groups 
are nicely centered around the corresponding transport sphere. The continued deleptonization after the neutrino burst 
caused a rather broad trough in the electron fraction profile. The region of neutrino emission coincides with the 
cooling region where the fluid is in weak equilibrium with the neutrino abundances. In the heating region, between 
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Fig. 4. — The last inelastic interactions of escaping neutrinos at 5 ms after bounce in the 13 Mq model. The abscissa of the graph is 
the radius, ranging from the center of the star to 80 km radius. The density markers at the bottom of the graph indicate the position 
of density decades, logjQ (p), where p is given in g/cm^. The thick solid line shows the electron fraction profile. The thick dashed line 
shows the electron fraction in weak equilibrium under otherwise unchanged conditions (same density, temperature, neutrino abundances, 
and spectra). These quantities belong to the ordinate at the left hand side. The total shaded area in the graph corresponds to the total 
electron neutrino luminosity at 5 ms after bounce in units of erg/s. The appropriate ordinate on the right hand side then carries the units 
erg/s/km. Note that the scale is 50 times larger than in Fig. Jsj. In this phase, almost all neutrinos escape directly after their production 
by electron capture on free protons (area below the white line). The emission of the neutrino burst occurs from a very narrow radius 
interval. This causes a steep drop in the electron fraction at the same position. The subdivision of the burst into different energy groups 
also shows a narrow energy spectrum of the emitted neutrinos. Inside 50 km radius, the trapped neutrinos are in weak equilibrium with 
the matter. In the burst region, the deleptonization time scale is only slightly slower than the shock propagation . Outside the shock front, 
it is much slower because of the low abundance of free protons. 



the cooling region and the shock front at 140 km radius, the reaction time scales are comparable to, or larger than, 
the infall time scale. Graph (b) shows the origin of the electron antineutrino luminosity. The electron antineutrinos 
are emitted fr'om a slightly smaller radius at densities exceeding 10^^ g/crn^ (while they were similar to 10^^ g/cnr" 
for the electron neutrinos). The electron antineutrinos decouple at smaller radii because of the higher neutron than 
proton abundance. The contribution of electron-scattered neutrinos (above the white line) is slightly larger than in 
the case of the electron neutrinos. The pair production process does not noticeably contribute to the luminosity at 
this stage. The equilibrium entropy (dashed line) is increasing with increasing radius. An infalling fluid element 
changes its entropy rather slowly by electron capture until it hits the accretion shock. At the shock front, most of its 
kinetic energy is converted into heat by shock dissipation. The heavy nuclei are dissociated, mainly into free nucleons 
which are good neutrino absorbers. However, before the stage at 50 ms after bounce, the entropy of an infalling fluid 
element has already reached or exceeded the equilibrium entropy, alone by shock dissipation. Only cooling is then 
possible during the continued infall. After tph ~ 50 ms, however, the equilibrium entropy is not reached by the shock 
dissipation and the fluid element indeed increases its entropy towards the equilibrium by neutrino absorption during 
its flight through the heating region (solid line in graph (b) 120 km radius). As the fluid element also tends to 
decrease the electron fraction (graph (a)), electron antineutrinos are preferentially absorbed. Once the equilibrium is 
reached, cooling becomes unavoidable if the fluid element is still infalling because the reaction rates at these densities 
are faster than the dynamical infall time. Convection in the heating region is expected to increase the neutrino heating 
efficiency (see e.g. the parameter study of Janka & Miiller (1996)), but not included in our simulations. We did not 
find any sign of shock revival in spherical symmetry with the included input physics. Graph (c) shows the origin of 
the /i- and r-luminosities. The /i- and r-neutrinos are produced at an average density slightly larger than 10^^ g/cm^. 
Almost no neutrinos escape directly from pair creation (area enclosed by white lines at the bottom of the figure). 
Most of the neutrinos have scattered off electrons before their escape (area enclosed by green lines). Pair production 
is not the dominant source of /x- and r-neutrinos. It has been shown, that the pro duction from bremss trahlung 
()ThomDson. Burrows, fc Horvath Il2000() and from electron flavor neutrino annihilation l|Buras et al. ll20'08a[) exceeds 
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Fig. 5. — The last inelastic interactions of escaping neutrinos at 50 ms after bounce in the 13 Mq model. The abscissa of the graph is 
the radius, ranging from the center of the star to 150 km radius. The density markers at the bottom of the graph indicate the position of 
density decades, logjQ (p), where p is given in g/cm^. The thick solid line in graph (a) shows the electron fraction profile. In graph (b) it 
is the entropy profile, and in graph (c) the electron chemical potential. The thick dashed line gives the equilibrium Ye in graph (a), the 
equilibrium entropy in graph (b), and the temperature profile in graph (c). We do not repeat the detailed explanation of the differently 
shaded and separated areas indicating the energy- and reaction-specific contribution of neutrino emissivities to the total luminosity of the 
star. Instead, we refer to the caption of Fig. JaJ or to the explanation given in the text. Graph (a) shows the origin of escaping electron 
neutrinos. Only electron capture contributes significantly (area below the white line). By definition, the region of neutrino emission 
coincides with the cooling region where the fluid is in weak equilibrium with the neutrino abundances. In the heating region, between the 
cooling region and the shock front at 140 km radius, the reaction time scales are comparable to, or larger than, the infall time scale. Graph 
(b) shows the origin of the electron antineutrino luminosity. The electron antineutrinos are emitted from a slightly smaller radius. The 
contribution of electron-scattered neutrinos (above the white line) is also slightly larger than for the electron neutrinos. By following the 
entropy profile from the right to the left (as an infalling fluid element would experience it) we flnd the expected abrupt entropy increase at 
the shock position. Behind the shock, the fluid element would drift more slowly inwards and neutrino absorption indeed leads to a small 
increase of the entropy in the region between 140 km and 100 km radius. However, the equilibrium entropy is declining towards smaller 
radii such that cooling becomes unavoidable once the fluid entropy has joined the equilibrium entropy in a still infalling state. Graph (c) 
shows the origin of the fi- and r-luminosities. Almost no neutrinos escape directly from pair creation (area enclosed by white lines at the 
bottom of the figure). Most of the neutrinos have scattered off electrons before their escape (area enclosed by green lines). 



the pair process production rate (both reactions are not included in our simulations). The latter reference finds that 
the /i- and r-neutrino luminosities show differences of 10% — 20% in the first 100 ms after bounce and converge to 
the standard luminosities afterwards. The spectra are not significantly different. This finding is also supported by 
our graph (c): The production site of the /x- and r-neutrinos is at a much smaller radius than their transport sphere. 
Hence, the neutrino luminosity is set by the (though energy-dependent) diffusivity between the location of neutrino 
production and the transport sphere. Moreover, graph (c) shows that the majority of escaping neutrinos scattered 
off electrons in their last reaction. Differences in the production-spectrum are likely to be washed out during the 
thermalization the neutrinos are experiencing while they are diffusing outwards to the transport sphere. 

Finally, we present the situation at 500 ms after bounce in Fig. (jSJ. This is after a long quasi-stationary phase 
of matter accretion and shock recession. The volume of neutrino emitting material has considerably shrunken with 
respect to the situation at 50 ms after bounce. But there are not much qualitative changes. The neutrinos with larger 
energies still escape from larger radii because of the corresponding staggering of the transport spheres in the steep 
density gradient at the surface of the protoneutron star. The shock has receded to a radius of 57 km. Graph (b) shows 
the origin of the electron antineutrino luminosity. The infalling fluid elements are now crossing the heating region 
that rapidly (with several thousand km/s) that there is no time for significant neutrino heating. This can be seen in 
the flat top of the entropy curve between 30 km and 50 km radius (solid line) . Even the cooling sets in with a slight 
delay and thermal balance is only reached at densities larger than 10^^ g/cm'^. Graph (c) shows the origin of the /x- 
and r-luminosities. Still very few neutrinos escape directly from pair creation (area enclosed by white lines at the 
bottom of the figure) . Most of the neutrinos have scattered off electrons before their escape (area enclosed by green 
lines). The continued cooling by /i- and r-neutrino emission becomes now clearly visible in the entropy profile at a 
radius of 18 km, where an entropy dip develops. At larger radii, however, between 25 km and the shock position, the 
temperature is rising and the electrons have become nondegenerate. The overlap of this material with the emission 
region of electron antineutrinos in graph (b) lets the emission of higher energy antineutrinos shift to lower densities 
than before. 

Figure ITJ shows the luminosities and rms energies of the neutrino flux after bounce on a longer time scale. In the 
13 Mq model, the luminosities decrease as a consequence of the declining accretion rate and continued deleptonization 
of the core. The electron flavor luminosities reach very similar values because the lifted electron degeneracy in a large 
part of the cooling region (see Fig. ((Sjac)) allows the electrons and positrons to be captured from similar chemical 
potentials. The luminosities are higher than the luminosities of the fi- and r-neutrinos because the latter do not have 
an accretion luminosity component. The rms energies show the usual hierarchy at the beginning, but after tpi, = 0.7 
s, the rms energy of the fi- and r-neutrinos falls below the rms energy of the electron antineutrino. This is also 
understood if one looks again at Fig. (|H|3c). While the emission of high energy electron antineutrinos is aided by 
shock-heated material settling at the base of the cooling region with moderate electron degeneracy, the layers around 
18 km radius, i.e. where the energy spectra of the fi- and r-neutrinos are set, are barely affected by the continued 
accretion on the still not very massive protoneutron star. This domain just slowly cools by neutrino emission (compare 
the entropy profiles in Fig. ^p) and (jHt))- The result are decreasing luminosities and rms energies of the fi- and 
r-neutrinos. The more massive 40 Mq model shows qualitatively different features. In order to understand them, we 
first discuss Fig. ©. Shown are the profiles of the mass flux through surfaces at constant radii in the two models. 
The dashed lines show the mass flux in the 13 Mq model. Nothing special happens there, the mass flux generally 
decreases in accordance with the decreasing density in the outer layers. Moreover, the contraction of the stiff core, 
which is far from its maximum mass, is minimal. A similar evolution is visible during the first 300 ms in the 40 Mq 
model. However, if we examine the massflux in the 300 ms time slice in Fig. ||2J) more closely, we find a variation 
in the density profile around 1000 km that falls in (from ^ 2000 km at 100 ms after bounce). The corresponding 
variation in the massflux or accretion rate leads to a step in the electron flavor neutrino luminosities between 300 
ms and 350 ms after bounce. The increase is of order 20%. The slope in their rms energies flattens slightly. More 
independent of the details of the progenitor model, however, might be that the protoneutron star in the 40 Mq model 
approaches its maximum mass much more rapidly because a high accretion rate is maintained when the outer layers 
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Fig. 6. — The last inelastic interactions of escaping neutrinos at 500 ms after bounce in the 13 Mq model. The abscissa of the graph is 
the radius, ranging from the center of the star to 80 km radius. The density markers at the bottom of the graph indicate the position of 
density decades, logj^Q (p), where p is given in g/cm^. The thick solid line in graph (a) shows the electron fraction profile. In graph (b) it 
is the entropy profile, and in graph (c) the electron chemical potential. The thick dashed line gives the equilibrium Yi^ in graph (a), the 
equilibrium entropy in graph (b), and the temperature profile in graph (c). We do not repeat the detailed explanation of the differently 
shaded and separated areas indicating the energy- and reaction-specific contribution of neutrino emissivities to the total luminosity of the 
star. Instead, we refer to the caption of Fig. \^ or to the explanation given in the text. Graph (a) shows the origin of escaping electron 
neutrinos. Only electron capture contributes significantly (area below the white line). Neutrinos with larger energies still escape from 
larger radii because of the corresponding staggering of the transport spheres at the top of the figure. The shock has receded to a radius of 
57 km and all regions are more compact. Graph (b) shows the origin of the electron antineutrino luminosity. The accreted fluid elements 
fall rapidly through the heating region without significant neutrino heating. Graph (c) shows the origin of the p- and r-luminosities. Still 
very few neutrinos escape directly from pair creation (area enclosed by white lines at the bottom of the figure). Most of the neutrinos have 
scattered off electrons before their escape (area enclosed by green lines). The continued cooling by p- and t- neutrino emission becomes 
now visible in the entropy profile at a radius of 18 km. At larger radii, however, between 25 km and the shock position, the temperature 
is rising and the electrons are nondegenerate. The overlap of this material with the emission region of electron antineutrinos in graph (b) 
enhances the emission of higher energy antineutrinos. 
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Fig. 7. — The luminosities and rms energies in the two models are shown on a longer time scale. In the 13 Mq model (dashed lines), 
the electron neutrino luminosity (thick line) and electron antineutrino luminosity (medium width line) converge. Both are much larger 
than the ^- and r-neutrino luminosities (thin line) because of the contribution from the accretion luminosity. The rms energies show the 
conventional hierarchy at the beginning. However, at 0.7 s after bounce, the rms energy of the ^- and r-neutrino luminosity falls below 
the rms energy of the electron antineutrinos. This is due to the fact that this low mass protoneutron star is quite incompressible with 
respect to the accumulated mass at the given accretion rate. The continued emission of ^- and r-neutrinos cools the deep layers (around 
r = 18 km in Fig. JHJ) independently from the dynamics of the outer layers. This is different in the 40 y\-odot mass model (solid lines). The 
protoneutron star comes closer to its maximum mass at 0.3 s after bounce such that it contracts appreciably with continued accretion. An 
increase in the accretion rates lets the electron flavor neutrinos step up. The accelerated transition from a mass accumulating stiff central 
object to a contracting compressible core affects the p- and r-neutrino properties. Shock-heated material is faster condensed to higher 
densities and the luminosities and rms energies of the p- and r-neutrinos start to rise continuously until the protoneutron star collapses to 
a black hole at ^ 0.5 s after bounce. 
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Radius [km] 

Fig. 8. — Mass flux through surfaces at constant radii in the 13 Mq model (dashed lines) and 40 Mq model (solid lines). The profiles 
are given at 0.1 s, 0.2 s, 0.3 s, 0.4 s, and 0.5 s after bounce. While the accreted mass piles up on the small protoneutron star in the 13 Mq 
model, it forces the massive protoneutron star in the 40 Mq model to hydrostatic gravitational contraction after tpi, = 0.3 s. An incoming 
variation in the accretion rate is visible outside of 1000 km radius. 



Table 1 

Enclosed mass at 100 km radius. 





13 Mq 


40 M0 


[s] 


[Mq] 


[Mq] 


0.0 


0.90 


0.88 


0.1 


1.19 


1.57 


0.2 


1.25 


1.77 


0.3 


1.28 


1.89 


0.4 


1.31 


2.02 


0.5 


1.33 


2.20 


1.0 


1.42 





"Shown is the enclosed mass at a radius of 100 km as a function of time for the 13 Mq and 40 Mq model. In the 13 Mq model, the 
accretion rate reduces quickly and the mass of the protoneutron star does not even come close to its maximum mass during the ~ 1 s time 
window so far explored with Boltzmann neutrino transport. The outer layers in the 40 Mq model are much denser. The accretion rate 
stays high and the simulation can be performed until the protoneutron star collapses. 



fall in. They have a significantly larger density in comparison to the 13 Mq model. The fast mass accumulation in 
the 40 Mq model becomes evident in Table (Q, which lists the enclosed mass at 100 km radius for different time 
slices in both models. As the accumulated mass in the 40 Mq protoneutron star gets closer to the maximum mass, 
the protoneutron star starts to contract faster by the general relativistic enhancement of the effective gravitational 
potential (an effect absent in Newtonian calculations). We observe in Fig. ((SJ that, after an initial decrease, the mass 
flux in the inner core is increasing again. This is by no means a "sudden" change on the short dynamical time scale 
of the protoneutron star. The contraction is a hydrostatic adaption to the accumulated mass in the gravitational 
potential. The change is, however, sudden on the time scale of the variations in the neutrino properties shown in Fig. 
JTJ. The /i- and r-neutrino luminosities and rms energies rise steeply. This happens because, by the contraction of the 
protoneutron star, electron-nondegenerate shock-heated material is condensed to densities where the main emission of 
heavy neutrinos occurs. We investigate the conditions at the locations of the main neutrino emission in more detail 
in Fig. ll^. We calculate the average conditions for the emission of a specific neutrino type according to Eq. ljB7|l . 
Instead of the average radius, we calculate here the average density (graph (a)), temperature (graph (b)), and electron 
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Fig. 9. — The thermodynamic conditions where escaping neutrinos have their last inelastic interaction are averaged with the weight of 
the contribution of the neutrinos to the total luminosity. Shown is the average neutrino emission density in graph (a), the average neutrino 
emission temperature in graph (b), and the average neutrino emission electron fraction in graph (c). Inelastic scattering of neutrinos off 
electrons is also included as an "emission" reaction. The typical conditions for electron neutrino emission are represented with a solid line, 
the typical conditions for electron antineutrino emission with a dashed line, and the typical conditions for fi- and r-neutrino emission with 
a dash-dotted line. The thin lines belong to the 13 Mq model and the thick lines to the 40 Mq model. The conditions are similar at 
bounce. They diverge afterwards. In the 13 Mq model, the main neutrino emission recedes to higher densities during the evolution of the 
simulation. In the 40 model, the main neutrino emission turns back to lower densities after gravitational contraction sets in at 300 
ms after bounce. The emission from the gravitationally compressed shock-heated matter at moderate densities competes more successfully 
with the emission from the cooler core at higher densities. 



fraction (graph (c)). These represent the typical conditions where an escaping neutrino makes its last energy-changing 
interaction with the matter. The conditions at the origin of electron neutrino emission are traced with a solid line, the 
conditions at the origin of electron antineutrino emission with a dashed line, and the conditions at the origin of fi- and 
r-neutrino emission with a dash-dotted line. We discuss first the 13 Mq model (thin lines). After an initial decrease 
up to the time of maximum neutrino heating around 100 ms after bounce, the average density at the sites of neutrino 
production increases steadily. This goes along with a temperature increase for all neutrino types. At least in the case 
of the very temperature-sensitive pair production rates for the /i- and r-neutrinos, however, the argumentation should 
be turned around: Because the temperature at a given density is slowly decreasing on the long time scale, the emission 
from lower density regions decreases and the average emission conditions shift to increasingly deeper layers, where 
the temperatures are higher. Graph (c) shows that the electron fraction at the place of emission is increasing for the 
electron flavor neutrinos. This is due to the rather high electron fraction in the shock-heated material behind the shock 
front, where the electron degeneracy is gradually lifted. The fi- and r-neutrinos escape from the floor in the electron 
fraction profile, which is steadily decreasing by continued deleptonization. This is reficcted in the declining electron 
fraction at the location of the production of heavy neutrinos in graph (c). The contraction of the protoneutron star 
in the 40 Mq model causes a qualitatively different characterization of the regions of main neutrino emission. The 
temperature increase by adiabatic compression is no longer fully balanced by neutrino cooling. The temperature at 
a given density starts to increase. The domain of the shock-heated material where the electrons are non-degenerate 
reaches down to deeper layers close to the place where the fi- and r-neutrinos are produced. The temperature increase 
makes these regions to significantly more efficient /x- and r-neutrino emitters such that the average emission density 
starts to decrease. In spite of this steep density decrease in graph (a), the temperature in graph (b) still increases at 
the average emission condition. The electron fraction in graph (c) rises dramatically as the mean neutrino emission 
region moves out of the Y,, trough in the electron fraction profile. As the production site of heavy neutrinos moves 
outwards to lower densities, the fraction of directly escaping neutrinos (without neutrino-electron scattering) also 
increases. Additionally, the forming "density cliff" shortens the escape path for few high energy fi- and r-neutrinos 
that may leave the star without thermalization. It can be expected that the not included reactions of fi- and r-fiavor 
neutrino production by bremsstrahlung and electron flavor neutrino annihilation will more significantly influence the /i- 
and r-neutrino spectra in this less opaque density regimes. Unfortunately, we cannot follow the evolution beyond the 
collapse of the protoneutron star because of the coordinate singularity at the formation of the Schwarzschild horizon. 
However, the neutrino luminosities are expected to decay on a short timescale (see e.g. l|Baumgarte et al. ..1996^') '). 

3. NUMERICAL IMPLEMENTATION 

In this section we motivate and describe our finite difference representation of spherically symmetric radiation 
hydrodynamics in numerical computations. The continuous equations of radiation hydrodynamics are unique. But 
the number of possible finite difference representations is large. We select ours according to the following requirements 
(with descending priority): 

1. In the limit of high resolution, the finite difference representation should converge to the continuous equations. 

2. The finite difference representation should lead to a numerically stable solution. 

3. The finite difference representation should not be less accurate than simpler schemes in regions where simpler 
schemes are sufficient. 

4. The computed evolution of the particle distribution function should imply an accurate evolution of expectation 
values like the particle density or the average particle energy. 

In our discrete ordinates (Sat) method, we discretize the full transport equation at once. This takes care of requirement 

(1) . In the first place, we choose the finite differencing such that the numerical scheme is stable according to requirement 

(2) . Moreover, we adjust it such that requirement (3) is met with respect to the important diffusion limit. In the finite 
differencing of the free streaming limit, we focus on accurate particle number luminosities and accurate particle energy 
luminosities. The accuracy of the angular moments of the radiation field is less essential in this domain because the 
neutrinos cannot contribute to the supernova mechanism in regions far outside the neutrinospheres. More accurate 
angular information in the free streaming li mit with the Sjy method w ould require a very large number of angular 
bins, an adaptive grid in angle space IfYamada. .Tanka. fc Suzuki 11199^ . or a postprocessing step with a ray-tracing 
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tool. All remaining degrees of freedom in the choice of the finite difference representation are used to optimize 
requirement (4). Of course, point (1) already guarantees that requirement (4) will successfully be approached if the 
resolution is sufficiently high. However, as we solve the transport equation at once, and with an implicit finite difference 
representation, resolution is computationally rather expensive. In the spirit of requirement (3), high resolution may 
only be required where the very details of the complete transport equation are relevant for the physical solution. In 
order to make low resolution results reliable, important physical laws must be represented accurately independent of 
the resolution setting. This includes the challenge of conserving lepton number and total energy in the simulations. 
Indeed, one might say that the power of a specific implementation of the Sat method is given by the extent to which 
requirements (3) and (4) can be satisfied at low resolution. 

A variable Eddington factor method (VEF) is a different numerical approach than the Sn method. There, re- 
quirements (2), (3), and (4) are satisfied by construction. A series of simpler radiation moment equations are solved 
and combined with an Eddington factor to produce the solution to the equations of radiation hydrodynamics. The 
quality of the variable Eddington factor method is determined by the capability to meet requirement (1). It de- 
pends on the interplay between the different moment equations and the accuracy of the Eddington factor (which 
may require the numeri cal solution of a model Boltzmann equation ( Burr ows et al. 1 120 00: Ra mpp fc Jankal '2002': 
iThompson. Burrows. Sz P into 2003,)). With careful implement ation, both methods should fulfill requirements (l)-(4) 
and therefore lead to similar results ijLiebendorfer et a'l~ll2003j) . 

3.1. Conservation Laws and Expectation Value Matching 

In the Sat approach, the complete dynamics in the phase space of the transported particles is determined by one 
kinetic equation (|15|1 . Its integration over the momentum phase space is equivalent to a continuity equation H24I) . its 
integration weighed with particle energies, e, is equivalent to an energy equation (|27|l . and its integration weighed with 
the particle direction cosine and energy would lead to a momentum equation. These derived equations certainly fulfill 
the macroscopic conservation laws as given in Eq. . But it is a challenge to obtain the same level of consistency in 
a finite difference representation of the Boltzmann equation. We illustrate this with the following little example: 

For simplicity we assume that a distribution function, f{t, x), of particles propagating with speed c only depends on 
time, i, and on one dimension in the phase space, x. The simplest form of a transport equation would then relate the 
time derivative of the distribution function to an advection term, 

-di^'d-x^^- ^^'^ 

If we now ask for the time evolution of the expectation value of the distribution function with respect to an operator, 
g(t, x), we apply integrations by parts to find 



(32) 



Here, the time evolution of the expectation value is described by other expectation values and the value of gf at the 
boundary of the domain of integration. This equation has an immediate physical interpretation: The change of the 
expectation value along the characteristics of the particle flow is simply given by the change of the weight along the 
characteristics of the particle flow because the phase space volume stays constant, 

+ fsf<"^- f(^. + '^)j<"^^ (33) 



cdt dx J J J \ cdt dx 

The second term on the left hand side is equivalent to the boundary term — [s/]^^, in Eq. (|32() . Let us now assume 
that numerical stability requires upwind differencing of the advection term in Eq. (|31|l . and, for simplicity, that the 
wind is unidirectional. A finite difference representation would then read like 

at dxi' 

We use the convention that a prime at the index points to the zone center value 1-1-1/2 while an integer index i 
points to the zone edge. We simplify the example once more and assume that the operator g is not time-dependent. 
With this, we track above integration by parts in the finite difference representation and calculate the evolution of the 
expectation value of operator g{x), 

9vh'dxe = g^' ^dx,' = ~ ^ g^' iU - .U-i) 

i—1 i—1 i—1 

n n— 1 



fi'dxi, - [gn'+ifn' - .go'+i'/o'] ■ (35) 
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' ^ dxe 
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Firstly, we remark that there is a discrete analogue to the integration by parts and that Eq. (|35|l has exactly the 
same structure as Eq. H32(l . Secondly, we note that the choice of finite difi^erencing for the advection term in Eq. H34|l 
determines the finite differencing of the weight function {gi'+i — Qi') / dxi' in the expectation value on the right hand 
side of Eq. (|35() . In the following, we will frequently use this type of investigation to optimize for requirements (3) and 
(4): First we choose a stable finite difference representation of a term in the transport equation. Then we evaluate 
relevant expectation values by performing discrete "integrations by parts" , analogously to this simple example. The 
emerging finite difference representations of the expectation values along one phase space dimension are then used as 
prefactors for the finite differencing of terms related to other phase space dimensions in order to build a consistent 
finite difference representation of the full transport equation. 

In order to construct a conservative finite difference representation of Eq. H15|) . we analyze the structure of the 
conservation equations in the continuum world, where insight is not buried in discretization indices. The emergence 
of lepton number conservation in Eq. (|24|l is straightforward, because the observer corrections are already written in 
a form that allows an immediate integration over energy or the angle cosine. Number conservation will also emerge 
naturally in the finite difference representation when we finite difference the same basic structure of the equation. 
The energy conservation equation in the frame of a distant observer, however, introduces the weight g{t^ a, ^, E) = 
E (T + ufi) in Eq. (|26|) . It depends on all four phase space dimensions and leads to the contribution of several 
expectation values of the distribution function in the evaluation of the energy conservation equation, 

J {r + u^i) E [Ct +Da + D^ + DE + OE + 0^- Cc] d^iE^dE = 0. (36) 

All terms in the Boltzmann equation (|15|1 . Ct, Da, D^, De, Oe, and O^, are written as the derivative of an expression 
with respect to time, rest mass, angle cosine, or particle energy. We now perform an integration by parts with 
respect to these integration variables x along the lines of Eq. H32|l . A multitude of correction terms proportional 
to d {r + ufi) E/dx arise. As in Eq. they describe the evolution of the weight q(t. x) = (T + uu) E alon g the 

characteristic of the particle flow. However, we showed in (Liebcnd orfer. Mezzacappa. fc Thielemamill2001|) that 
(r + uu) E is nearly a constant of motion along the characteristics. Therefore, it is natural that most of the partial 
terms d (T + ufi) E/dx actually cancel in the total evolution of the radiation energy in the frame of a distant observer. 
If, like in Eq. H35|l . the integrations by parts are also carefully followed in the finite difference representation, mutual 
cancellations of important expectation values can be forced to be exact — independently of the resolution. Because 
these canceling terms individually reach large values around and after bounce, this expectation value matching is an 
essential step for the conservation of energy in the finite difference representation. We perform the integrations by 
parts in Eq. H36|) and start with the identification of canceling contributions in the following overview of contributions 
from Ct, Da, D^, De, Oe, O^, Cc in the transport equation (|15|l : 
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We will neglect terms that are nonlinear in the radiation field, i.e. t erms that describe a gravitational interactio n of 
the radiation field with itself. From the hydrodynamics equations in ijMisner fc Sharp1ll964l: IMav fc White ill 966(1 we 
then derive the following useful relationships for the coefficients in Eq. H37|l . 
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We may now identify cancellations in Eq. H37() and label them for later reference: The fourth term of Da cancels with 
the first term of Ob {D^O^^) by Eq. H39(l . The first and second term of cancel with the third and fourth term 
in Ob (Dj^O^). These are the only 0{v/c) cancellations and therefore the most important ones. In an analogous 
notation, we find the following higher order cancellations: {CfD^) by Eq. {D^D^), {0%0l), (Of OjJ"*). The 

remaining terms, {CfD^D^Oj^), reduce by Eq. H4U|) and the definition of F to the general relativistic term. 

They enter in this form the radiation energy conservation equation H27|l . 

We will adopt the following strategy in the finite differencing of Eq. H15|l (the reader is invited to draw lines in Eq. 1)3 7|l 
to visualize the relationship between the different terms): (i) The finite differ encing of Ct and are straightforward, 
(ii) Appreciable experience has been gathered in previous work (see e.g. llLewis fc Miller I '1984')') with the finite 
differencing of the 0{v/c) terms in Da and D^. They are assumed to be well-chosen in |Mczzacappa & BruemJ 
Il998aj) and therefore not subject to changes, (iii) Based on this, the cancellation (-D^O^) dictates the finite difference 
representation of {d\n p/adt + 2u/r) in 0|; and O^. The cancellation (O^O^) sets the finite difference representation 
of {din p/adt + 2u/r) in O^ and O^. (iv) The cancellation {Djj^O^) dictates the finite difference representation of 
u/r in Ob* and Of. The cancellation (Of OjJ"*) propagates the finite difference representation to u/r in O^^. (v) The 
evaluation of {Cf DaD^Oj^) according to Eq. 1)41(1 constrains the finite difference representation of Td^/dr in D^ and 
therewith defines D^. (vi) And finally we choose a different finite difference representation for Tud^/dr in D^ and 
I?^, as suggested by the cancellation (CfD^'^D^). By these six chains, all terms in Eq. (|15f) become constrained. 

3.2. Adaptive Grid 

AGILE-BOLTZTRAN has the Capability to dynamically adapt the computational grid to the evolution of the stellar pro- 
file. Each radial zone contributes to the computationally expensive implicit solution of the Boltzmann equation. The 
adaptive grid allows the concentration of these zones to domains where radiation transport is actually active and impor- 
tant. It also allows a pro per shock capturing and propagation . The fundamental ide as of the dynamically adaptive grid 
have been developed in ijWinkler. Norman, fc Mihalas 11198^ iDorfi fc Drurv Ill987t) . The first reference describes how 
the continuous equations of radiation hydrodynamics are formulated on an adaptive grid and the second reference pro- 
vides a numerically stable prescription how to adapt the local grid point concentration to steep gradients in the solution 
vecto r. The details of our implementation of the adaptive grid are described in (Liebendorfer, Rosswog & Thielemann] 
l2f)02|l . Here we resume the key points and describe the extension of the approach from the hydrodynamics code agile 
to the radiation quantities. 

We concentrate on observers at rest in their slice. They reside in local orthogonal reference frames which we assume 
to be collinear with global coordinates A. Consequently, coordinates A have vanishing shift vectors everywhere. We 
introduce another system of global coordinates, B, having arbitrary shift vectors on the same time slices. Let us select 
fixed grid labels in coordinates B: {qi g R, z = 1 . . .n}. The fixed labels in coordinates B define a moving grid in 
coordinate system A: {a{t, Qi)}- The grid labels in coordinate system B cut the continuum into zones Aq on each time 
slice. These zones contain a time dependent selection of observers in coordinate system A: Aa ~ {a{t,q),q e Ag}. 
We define a zone integral of the observations 

{y{t,a)) = / y{t,a)da (42) 

JAa 

that is based on the observations y of single observers a in the zone. In coordinate system A, the borders of the zones 
change with time according to the shift vectors in coordinate system B. In Newtonian parlance, one can identify 
the moving zones in coordinate system A with the motion of an adaptive gr id. The important point is th at the 
adaptive grid only regroups the observers into new zones. As stated already by l|Winkler. Norman, fc Mihalas i ri9841. 
the adaptive grid never changes the reference frame of the observations. It solely determines the width and location 
of a zone on the time slice for the zone-integration of the observations y in the comoving frame A. 

In order to formulate time evolution equations, we are interested in the time change of the zone-integrated quantity 
{y{t, a)) in a specific zone Aq. Taking into account that the boundaries of the integration over Aa are time-dependent, 
this leads to 

d , dy{t,a) , 

— / 2/da= / — 17 — da + 

The notation d/dt is used for time derivatives at constant grid label q while Lagrangian derivatives, d/dt^ are evaluated 
at constant baryon mass a. The left hand term is the time change of the zone-integrated variable y between fixed grid 
labels in the numerically accessible coordinate system B. The first term on the right hand side is the time derivative 
in orthogonal coordinates A as used in the physical equations, e.g. on the left hand side of Eqs. ||HJl-lIHl)- The second 
term on the right hand side of Eq. (|43|l is due to the motion of the grid. It is evaluated at the zone boundaries and 
accounts for the observers that fall into or drop out of the zone during an infinitesimal time step. The rate of coordinate 
change at a fixed grid label, u"'^^ — ~da{t,q)/dt, defines a "grid velocity" with respect to coordinates A. Note that 



da{t, q) 



dt 



(43) 
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the adaptive grid corrections in Eq. H43() enter in a conservative form of fluxes from zone to zone. If the time change 
of the variable y is integrated over the whole domain, all grid corrections at intermediate zone boundaries cancel. The 
adaptive grid does not affect conservation properties of the variable y. This important feature is maintained in the 
finite difference representation of Eq. (|43() , 



At 



dy\ Oi+i - flj+i a, 
-lAa,+y.^, V^' 



At 



(44) 



The overbars in this equation mark the quantities that are evaluated at the old time, before the time step At. 

The application of Eq. H44|) for the evaluation of Lagrangian time derivatives on the adaptive grid is usually 
straightforward. However, a side-remark concerning our ad hoc "burning" of silicon to nuclear statistical equilibrium 
(NSE) at the edge of the iron core should be made here. In a Lagrangian scheme, we would simply convert zone i + 1 
to NSE during infall as soon as the burning time scale becomes considerably smaller than the evolution time scale. 
The conversion from a given composition to NSE is handled such that energy remains conserved with the respective 
equations of state. If the adaptive grid is applied to the energy equation, it automatically performs this conversion 
also to the continuous mass advection between zone i and i + 1. This is perfect if the grid moves outwards such that 
advected silicon is continuously converted to NSE. It is less perfect if the grid moves inwards, because this leads to the 
advection of material labelled with "NSE" from zone i to "silicon" in zone i + 1 until the whole zone i + 1 carries the 
label "NSE". At this point the whole zone is converted back to NSE. Without additional measures, the advected NSE 
material would therefore transiently be treated as "silicon" . By energy conservation, the temperature of zone i + 1 
would drop because energy would be invested to unphysically "unburn" the advected NSE material. We eliminate this 
effect by externally feeding the zone with the difference in the nuclear binding energy between advected and internal 
material such that no sudden changes in the thermodynamical state will take place. As soon as the zone is converted 
back to NSE, we extract the same amount of energy from the burning energy to restore the balance. 

If the adaptive grid technique is applied to radiation quantities that depend on momentum phase space variables 
like the particle energy, E, or the angle cosine, fj,, further considerations are required. A particle, that is observed in 
the comoving frame of zone one with an energy E and a propagation angle fj,, is observed in the comoving frame of an 
adjacent zone two with energy E' ^ E and ji' ^ ^jl because of frequency shift and angular aberration. If the moving 
zone boundary sweeps over the particle such that it is assigned to the adjacent zone two, the momentum phase space 
state of the particle in zone two should be changed to E' and in order to describe the same physical state the particle 
had when it was a member of zone one. The motion of the zone boundary has a second effect: The location of the 
center of zone two moves towards zone one and fluid momentum is advected from zone one to zone two. This means 
that the particles that do not change the zone are as well subject to frequency shift and angular aberration corrections 
in order to maintain their exact original state after the redefinition of the zones. In the more general derivation of the 
adaptive grid equations in Winkler and Norman (Winkler. Norman. & Mihalas 1984), these two corrections do not 
appear. In order to convince ourselves that they cancel, we calculate them explicitly. 

We assume that we have a global description of the particle state in the momentum phase space that does not 
depend on the location of the particle or the velocity of the fluid element. We call these global parameters "impact 
parameter", 6, and "energy at infinity", e, because the quantities 



b = 



r + M/i 

e={T + uii)E. (45) 

are in good approximation constants of motion along the general relativistic particle trajectory 
l)Liebendorfer. Mezzacaopa. & Thielemann! I2001D . Once the description of the particle momentum phase space 
is invariant, it is safe to reassign neutrinos between zones according to the adaptive grid equation We substitute 
the general variable y{t, a) by the distribution function _F(i, a, 5, e). 
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Evaluating the first term by the chain rule gives 

d_ 

dt 

fa(92) rgpdn 
a(gi) [dfi dt 
The second term can be written as 



1(92) 

a(9l) 

dF dE 
dE'dt 



■da 



Fda= — 



1(92) 
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Finally, the term on the right hand side evaluates to 



a(qi) dt 



—J— da - 

a(qi) dt 



£1(92) 



a{qi) 



dF dn 
dfi dt 



dF dE 
dE~dt 



da. 



(49) 



Substitution of expressions (|47|l - l|49(l into Eq. H46|l shows that the correction terms indeed cancel. Hence, we can as 
well directly apply Eq. (|43|l to the neutrino distribution functions, F {t, a, /x, E), defined in the comoving frame 



d_ 

dt 



0(92) ra.{q2) 

Fda ~ / —:-da 



a(qi) 



a[qi) 



dt 



F 



,9a 
di 



(50) 



We note however, that energy conservation in the adaptive grid corrections will become resolution dependent if we finite 
difference Eq. 1)51) |l instead of the much more complicated Eq. H46|l . Particle number conservation is still guaranteed 
at machine precision. We quantify this effect in our energy conservation test in section [4.61 

The implementation of an adaptive grid consists of two parts: We discussed above how the physical equations are 
evolved on a dynamically moving grid. Here we summarize how the motion of the grid itself is controlled. On the 
one hand, it has to increase the local grid point concentration where steep gradients in the solution vector need to be 
resolved. On the other hand, it has to capture self-similar flows and to propagate them through the computational 
domain. If there are traveling features in the flow, it is advantageous if the adaptive grid can just change the 
zone locations instead of significa ntly changing all the physical variables describing the state in the zones. A useful 
grid equation has been found in ijDorfi fc Drurv lll987|) . We transferred the approach from Eulerian coordinates in 
Newtonian space to comoving coordinates in general relativistic time slices such that we retrieve a Lagrangian scheme 
if the adaptive grid is switched off. The grid equation of Dorfi and Drury is based on a resolution function, 



R{t,a) = [l + Y. 



0"^=' dy{t, a)" 



y 



da 



(51) 



where the sum includes a selection of relevant variables y. Currently, our grid position is guided by gradients in the 
velocity and logarithmic density profiles. The velocity contribution focuses on shocks and the density contribution 
draws grid points to steep density gradients where the resolution is important for the radiation transport. We set 



the global scale in this equation according to the total variation, y' 



,scl 



concentration, n 



-,scl 



If we define a grid point 



to the requested resolution l|5T]l. 



After the substitution of 



{da{t, q)/dq) , the grid equation requires that the the grid point concentration is proportional 



d_ 

dq \R 



R 




(52) 



(53) 



we note that this presc ription enforces a constant generalized arc length per grid point interval Ag. In the same way as 
iDorfi fc Drurv I 1)198 7^. we apply operators for space smoothing and time retardation to the grid point concentration 
before solving Eq. H52|) implicitly with the hydrodynamics equations. As the adaptive grid requires corrections in the 
evaluation of time derivatives, we have payed attention to reduce the occurrence of time derivatives to a minimum 
in the implementation of the Boltzmann equation on the adaptive grid. Therefore, we can relax the notation in the 
following: We will use the notation d/dt for the time derivatives at constant enclosed mass. This makes the discussion 
of the radiation transport equation more consistent because the enclosed mass is just one of four dimensions in the 
phase space, and the derivatives will always be partial. The few time differences at constant zone labels will explicitly 
be written out in the finite difference expressions. 

3.3. Finite differencing of the transport equation 

We split the description of the finite differencing of the Boltzmann equation H15|l into separate terms according to 
the list in Eq. (|15|I - H22|) . Each term is prepared in its own subsection. In section we describe how all the individual 
terms are gathered to form one implicit solution of the complete transport equation and an operator split implicit 
solution of the hydrodynamics equations. 

The knowledge of a complete set of primitive variables at a certain time allows the derivation of all other quantities of 
interest. Our choice of primitive variables is listed in table The discretization in space is indicated by the index i. 
We adopt the convention that integer values of i point to zone edges while i + 1/2 refers to the mass center of the zone 
between edge i and i + 1. Wc abbreviate half-valued indices with a prime, i' = i + 1/2, in order to not further challenge 
the readability of the finite difference equations. Of order 100 spatial zones are used on a dynamically adaptive grid. 
The particle distribution functions, -Fi'j'.fc', additionally depend on the particle momentum phase space, i.e. the 
angle cosine of the particle propagation direction, /i, and the particle energy, E, measured in the comoving frame. 
Both momentum phase space dimensions are discretized with zone center values, /ij' and Ek' , in compliance with the 
above-mentioned conventions. We use Gaussian quadrature in the angles fij' with weights Wj>. They are derived from 
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Table 2 
Primitive variables. 





primitive variable " 


primitive variable 




ai 


enclosed rest mass 


Pi' 


rest mass density 


ri 


radius 


T r 


temperature 


Ui 


velocity 


Y -t 

e,z 


electron fraction 


nii 


enclosed gravitational mass 




lapse function 






Pi',j',k' 


neutrino distribution function 



"The primitive variables are listed with an index i if they live on zone edges and an index i' if they live on zone centers. The momentum 
phase space for the neutrino distribution function is labelled with an index j' for the angle cosine of the propagation direction and an index 
k' for the neutrino energy. 

Legendre polynomials and normalized such that j\dfi = Wji = 2. The range from inwards to outwards propagation 
is currently resolved with 6 different propagation angles, but we are free to choose a larger number of directions for 
higher resolutio n (see section [4. 5|l . The energy grid is set by geometrically increasing zone edge values Ek- Following 
iBruenn 1 1)2002|) . we define the zone center energies by 

I ^k+i + ^k+iEk + El 

Ek' - y ^ . 

This choice implements the phase space volume El,dEk' = (^'l+i ~ E^^ /3 in an exact manner. The zone center 
energies also form a geometric series. We currently use 12 zone center energies ranging from 3 MeV to 300 MeV, 
but this number can also be increased for higher resolution (see section 14.511 . Thermodynamical and compositional 
quantities are obtained from the primitive variables p, T, and Y,, by the equation of state. Properties of the particle 
radiation field are obtained by the evaluation of the expectation value of the corresponding particle property with the 
primitive particle distribution function F. 

3.3.1. Time derivative of the particle distribution function 
The time derivative of the particle distribution function, Eq. Hlt)|) . is implemented in the following way: 

r - ^ f Fi,^j,^k'dai' ~ Fj^^y^k'daj, ^ ^ „roip* l\ ftrA\ 

The quantities with an overbar refer to the previous value of the variables. All other quantities refer to the updated 
values after the time step dt ^ t — i. On a Lagrangian grid, the mass in a spherical shell, da,' — a^+i — a^, is a constant 
and the relative velocity (in mass) between the matter and the grid = — (a^ — a^) /dt vanishes. In this case the 
expression would simply reduce to Ct = a^^ {Fi'j'^k' — ^i' j',fc') /dt. On the adaptive grid, however, we apply Eq. 
and obtain Eq. (|54ll . The additional particle fluxes across moving zone edges, uj'^^'^F*^-, j., , are evaluated by upwind 
differencing in order to guarantee numerical stability: 

F*; , ^ / ^''-iJ'.fc' if w^' > 
hj',k' I Fiij'^k' otherwise. ^ ' 

Note that this particle advection is determined by the grid velocity alone, and not by any physical particle speed. The 
grid velocity can in principle become arbitrarily large. Unlike the advection due to particle propagation, it does not 
depend on the particle mean free path. 

3.3.2. Particle propagation in space 

The term Da in Eq. (|17|l accounts for the particles that are propagating into and out of a spherical mass shell. 
From first considerations it is evident that this term has to be proportional to the dot product between the particle 
propagation direction and the shell boundary, i.e. the angle cosine /i. By this property, the contribution becomes 
proportional to the gradient of the first angular moment of the particle distribution function when the Boltzniann 
equation (|15|l is integrated over the propagation angles. In the free streaming limit, the advected particle number in a 
time step (always in the perspective of a comoving observer) is large with respect to the particle number in the zone. 
Upwind differencing of the advection terms is appropriate to limit destabilizing errors in the advection fluxes. For 
discrete angle cosines, /ij' , the direction of the particle "wind" is simply given by the sign of fiji (A particle with fj, — 
propagates tangential to the mass shell). In diffusive conditions, an asymmetric differencing can lead to an overestimate 
of the first angular moment because of improper cancellations among the contributions of the isotropic component 
of the particle field. Based on transport coefficients Pi^k' , we therefore interpolate betwee n upwind differencing in 
free streaming regimes {l3i,k' = 1) and centered differencing {Pi^k' = 1/2) in diffusive regimes ijMezzacappa &: Bruenn I 
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Il993aj) . The chosen transport coefficients, however, have the disadvantage, that they are not perfectly one half in the 
diffusive regime inside the neutrino sphere. This can cause a spurious numerical contribution to the diffusion flux in 
the hydrodynamic limit where the particle density, multiplied by the propagation speed, exceeds the particle flux by 
orders of magnitude. In order to satisfy the equilibrium constraint (73) in (Mezza cappa fc Bruenn 1993aj) everywhere 
inside the neutrinosphere, we use transport coefhcients according to 

p^^,^f 1/2 _^ if 2dn>X,, 

\ {2dri/ Xi^k' + ^) otherwise. 
That is, the transport coefficients at zone edges depend on the ratio between the distance dvi between zone cen- 
ters and the neutrino energy dependent mean free path A,;.fc' . We finite difference the propagation term as in 
(fMezzacappa fc Bruenn ■■199 3al : 

Da = — [4:nrf_^_j^ai+ipi+iFi+ij'^k' - ^Trrfa^piFi^j, ^k'] (57) 

with 

C^iPiFij'^k' — A,i',fc'aj'-lPi'-l-Fi'-l,j',fc' + (1 ^ Pi,j',k') Cti' Pi' Fi' J' ^k' 

for outwards propagating particles {fij' > 0), and 

c^iPiFij'^k' = (1 ^ /3i,i',fc') Q^i'-iPi'-i^i'-ij',fc' + Pi,j',k'Cti'Pi'Fi'j'^k' (58) 
for inwards propagating particles {pj' < 0). 

Based on this finite difference representation, we calculate now the expectation value D'^ in Eq. H37|) for a later 
match with the observer corrections described in subsections 13 . 3 . 4l and 13 . 3 . 5l For we identify in Eq. H37I) the weight 
Ui+iiij'Ek' in front of the space derivative in Eq. l(57|) . As in Eq. we perform an integration by parts over da in 

the finite difference representation, 

Mi'^fe' r 2 2 1 

|.47rrj_|_iQ;i+iPi+iUi+iFi+ij'^fc' - Anr^a^piUi+iFi^j' ^k' \ 



ai'dtti' 

ai'daii 

Airrf Ui+i -Ui 2 

-Hj'Fk'aiPiI'ij'^k 



ai' dui 

If we would additionally integrate over the momentum phase space, with E'^, as the usual measure of integration, 
the last term would correspond to = —Airr^p^K. We insert in this term the interpolation prescription H58|l for 
ctiPiFiji ^k' on zone edges and collect the two contributions with the zone center index i' in the distribution function 
to obtain a finite difference representation for A = —Anr'^p/T ■ du/da: 

Ai',k' = ^7 {"^i+l (■"i-l-2 - Ui+i) (1 - f3i+l,k') + rf {Ui+i - Mj) A,/c') 

if < 0, and 

A-t'.k' = r= — (r^+i (wj+2 - Ui+i) /3i+i,fe' + r,^ (mj+i - Ui) (1 - A.fe')) (59) 

if fij' > 0. We use exactly the same analysis to evaluate the finite difference representation of = — 47rr^p^iJ in 
Eq. (|37|l . The weight for this term in front of Eq. H57|l is Ti+iEk'. The result is a finite difference representation for 
A = -Anr^p/T ■ dV/da: 

At: p. 

iW,k' = J 

if pj' < 0, and 



Aj',fc' = :fT— ^7 — (rf+i (r,:+2 - Ti+i) (1 - (3i+i,k') + rf (E^+i - r,)/3,,fc') 



A.^fc- = -F^^^ {r-f+i (r«+2 - r,+i) p,+,,k' + rf (r,+i - r,) (i - ^,^0) (60) 

otherwise. These definitions allow to write the contribution of space propagation to the total energy conservation, 
J2j,k (^i+i + Ui+ipf) Da.i'W-j'El.dEk', in the form of Eq. 

D^'^^ : — ^ — [A-Krlj^^ai+ipi+iTi+iHi+i - AurlaiPiTiHi] 

OLi' CL(Xi' 

+ ^ i^i',k'Fi'j'^k'Pj'Wj' EI, dEk' 

+ Y,^i+iAi',k'Fi'j>,k'PfWj,EldEk'. (61) 
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3.3.3. Angular advection from spatial propagation 

If a particle propagates in outwards direction, its direction cosine increases towards unity. Hence, particle propagation 
causes angular advection relative to a fixed grid of particle direction cosines fij' . This contribution to the time evolution 
of the particle distribution function is described in the term Z)^ in Eq. (|18() . The finite difference representation of this 
term is chosen such that an equilibrium between an iso tropic radiation field and stationary matter remains undisturbed 
(jLewis fc Miller 11X9841 iMezzacappa fc Bruenn iri993aj) . Following this approach, we set 

- (rm !^['r~1j - G'^i ) ^ (0+i^^.',,+i.fc' - QF,,,^k') . (62) 
V 2 - rf\ J 

Characteristic for the general relativistic equation is the gravitational term ~ Td^/dr in Eq. (|62|l . It implements 
gravitational bending. We outline its finite difference representation below in Eq. (|68|) . The differencing of the 
coefficients, C, = 1 — ji^ , is defined by 

0+1 - = -2^j'Wj'. (63) 

The angular integration of the term Z?^ produces the zcroth and second angular moments of the particle distribution 
function. Its finite difference representation is therefore not as sensitive to cancellations in the diffusive limit as in the 
case of the spatial advection term Da- Upwind differencing can be justified. The angular "wind" always points towards 
^ = 1. However, for reasons of completeness and consistency, we also use centered differencing in the diffusive regime. 
With angular transport coefficients 7i',fc' := Pi',k', we interpolate the values of the neutrino distribution function on 
angular zone edges by 

Fi'.jM' = li'.k'Fi'j'-i^k' + (1 - li\k') Fi- ^k' ■ (64) 

Again, the terms D^^'^ in Eq. (|37|l arise from an integration by parts of the Boltzmann equation (|15|l with weight 
iXi+i + Ui+iUji) Ek/ ■ But this time the integration is over the angle cosines, /i. The part with the weight Ti^iEk' 
does not depend on the angle cosine, its contribution vanishes. For the calculation of the other part we introduce the 
abbreviation 

for the prefactor in Eq. Ht)2|l and perform the integration by parts again in its finite difference representation: 



3,k 



u. 



iCj+it^j' Fi' j+i^k' — CjtJ-j'Fi'j,k')Ei,,dEk' 



— ^^Uj+iTTi+i {Q+ifiji Fi' j^i^k' — CjlJ-j'-iFi'j^k') E^idEk' 

— ^ Ui+iTi+i(jFi> j^k' (a^j' ~ Mj'-i) EydEk' 
= — ^ Ui^iTi^iFiij>,k'El,dEk' 

j-k 

X [le.k'Cj+i (Mj'+i - Aij') + (1 - li',k') Cj (Pf - Mj'-i)] ■ (65) 
The expression in the second line vanishes because of mutual cancellations in the sum and Qi — Cjmax+i ~ 0- With 
C = 1 — /i^, we can identify the result as a special finite difference representation of the terms D^^ and D^^ in Eq. 
1)37(1 . Following the strategy of the previous paragraph, we extract from Eq. ((65(1 a finite difference expression for the 
expression B = (l — u/r, 



Bi',j',k' = 7;-T^ [li\k<]+i (Aij'+i - Mj') + (1 - li'.k') Cj (Mj' - Aij'-i)] • (66) 



With this definition, we write the contribution (Fj+i + Ui+ijij') D^^i'Wj' E^,dEk' to the total energy evolution in 
the form of its continuous analogue in Eq. 1(3 7() . 

: - r,+i ^ B,,,j,,k'F,,,j,,k'Wj,El,dEk' + u,+iGf+i ^ F,, ,k' E^dEk' 

j-k j,k 

X hi',k<3+i (Aij'+i - t'-j') + (1 - 7i',fc') ilH' - Aij'-i)] • (67) 
The finite difference representation of G^^^ for the gravitational bending still needs to be defined. We make its 
definition depend on the finite differencing of = Td^/dr in Eq. 1(84(1 . which will account for the gravitational 

)|; in Eq. (|?7|) it contributes with m+iGf^^ 



frequency shift of the particle energies. In the term in Eq. ((37() it contributes with Ui+iGf^j^Kii to total energy 



conservation. If we define 



^'+i-2^r.+i-r. -2^'+^' ^^^^ 
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this will support the desired cancellation Cj + Dj^ + = under two conditions: (i) The hydrodynamics equations 
should accurately implement Eq. (|38|l . i.e. 



If this is not the case, an alternative discretization based on the time derivative of F as a function of the velocity 
and the gravitational potential has also provided satisfactory results (see Eq. for the evaluation of du/dt on the 
adaptive grid), 

(ii) The particle distribution function should be close to isotropic. In the isotropic limit we can use the relation J = ZK 
and obtain from Eq. H65() a simplified contribution to energy conservation 



k j 



Note that the Gaussian quadrature in the angle discretization provides the exact integral for J fj^d^ = 2/3. Both 
conditions are fulfilled in the high density domain of the protoneutron star where discretization errors most dramatically 
affect energy conservation because of the large radiation energy densities. Of course, Eq. H68(l provides also a valid 
finite differencing of gravitational bending outside of the protoneutron star, where the effect is physically important. 
But the finite differencing is not specially tuned to enforce the cancellation C"^ + D^^ + = outside of the diffusive 
regime. 

3.3.4. Frequency shift from observer motion 

We continue our description with the observer correction in Eq. (|20|l . We start with a review of the finite 
differencing procedure presented in l|MezzacaDDa fc Bruenn1ll993aj) in order to extend it from local to global energy 
conservation. Global energy conservation requires the cancellation of observer correction terms in Eq. H37|l with terms 
from space propagation and angular advection. Therefore, it is convenient to express the observer corrections in terms 
of the same physical quantities A and B we have used before, 

^ f d\\ip ^ 2u^ 
\ adt 

^ = (1-/^^)7- 

Oe = {,'A-B)±^[E^F]. (70) 

In the following, we solve the equation for the frequency shift correction, Ct + = 0, by the method of characteristics. 
We convert the time derivative at constant energy to a time derivative along the characteristic in the energy dimension of 
the momentum phase space. This eliminates the partial derivative with respect to energy and facilitates a conservative 
finite differencing of O^;. Note that the characteristic used in this section implements the frequency shift c orrection 
alone, it differs from the free propagation characteristic prescribed by the full Boltzmann equation. As in ijBruenn I 
^85; Mczzaca ppa fc Brucnn 1993a.) we write the prefactor of the correction as a time derivative of a quantity Rf = 



aini? 



adt 



such that Ct + Oe — becomes 



It is now possible to transform from the "Eulerian" variable, x = E, to a "Lagrangian" variable, y = E/ Rf, along the 
characteristic by the chain rule (d/dt)^ ~ {d/dt)y + (dy/dt)^ {d/dy) (the subscript of the bracket denotes the variable 
that is kept constant for the differentiation). Eq. (|71|l simplifies to 

d r„,^n\ dRf ^ „ 9 



E/Rj 



^ rpSpA f d[E/Rf] \ d[E^F] f d 3 
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For a small section of the energy phase space E'^AE = (£'| — Ef) /3, this relationship leads to 

d 



dt 



[e'^fae] 



= 0. 



(72) 



E/B,f 



The validity of Eq. H72|) for arbitrary distribution functions F leads to the following interpretation: The observer 
correction shifts the particles that initially reside in the energy interval E^AE along the characteristic with constant 
E/Rf in the energy phase space. This allows us to determine the evolution of any other particle property in analogy 
to Eq. (|33|l . For example, the specific energy of the particles in this phase space interval, de ~ E^FAE, evolves 
according to 



[E^FAE] 



E/Rf 



E/Rf 



dlnRf 
dt 



de. 



(73) 



A finite difference representation of Eqs. H72|l and H73|l has been given in ijMezzacaDpa fc Bruenn lfl993aD . Consider a 
particle energy group k' , with a neighbor group k' + dk, dk = ±1. Eq. H72|l tells us that the number of particles before 
the correction. Fit jr ^k' E'^idEk' , is equal to the number of particles after the correction. The distribution after the 



correction is represented by a diminished number of particles Fi' ji ^k' E"^' dE^' 
number of particles n^^ ji k'+dk neighbor group k' + dfc. 



n^i ■, J., in group k' and an additional 



F, 



EtdEk 



F 



i',j'.k'J^k 



:'Ei,dEu 



+ 



n , 



i' ,j' ,k' -\-dk 



(74) 



Eq. H73|l now defines a similar relationship for the particle energies 



Fi'j'k'EuidEk 



(F,,,,,,k'El,dEk' ~Ek.n^,^^,^k) 

- {^l]>A^,^k' ~ B,,j,^k') F,,j,^k'EldEk'a^'dt, 



(75) 



where Ai'^k' and Bi'j>^k' stand for a finite difference representation of Eq. H70|l . Eqs. (I74() and (|75|l uniquely define 
the solution 



i' ,j',k 



dEk' 



E, 



k'+dk 



-E, 



-El,F,,,j>^k'ar>dt 



''^t',j',k' —^i' ,j' .k'-dk^ 



(76) 



which leads, by the update Fiiji^k' = Fii^i^k' + {j^t' j' k' ~ ''^e j' k'^ / {E^idEk'), to the following finite difference 
representation of the frequency shift term in the Boltzmann equation H15|l : 

1 



El dEk' 

t^^j'Ai'^k'-dk — Bi/ J' ^k' -dk) 



dEk' 



k'-dk 



Ek' — Ek: 



k' — dk 



-Ek'-dkEt',j',k'-dk 



i' ,k' 



Bi',j',k') 



dEk' 



-EuiFi 



k'^i' ,i' ,k' 



Ek'+dk — Ek' 

Finally, we calculate the contribution of the frequency shift to energy conservation. 



(77) 



(r^+i + Ui+i^ij^) OE,i'Wj' El,dEk' 



The summation over k in Eq. H77|) with weight Ek' and measure of integration E^ can be simplified with a discrete 
"integration by parts" as in Eq. H35|) . If we neglect for the time being the boundary terms because of a small Ey in 
the lowest energy group and a small Fi'^ji^k'^^^ in the highest energy group (a correction for these terms will be made 
later in subsection 13. 3. 7|l . the energy contribution from Oe becomes 



O 



1-4 



— ^ (Ti+i + Ui+ifiji) (fi'^/Ai'^k' — Bi'j'^k') Fi>j>^k'Wj'El,dEk'. 



(78) 



If we use the finite difference representations (|59|) and H66|l for Aii ^k' and Bi' ji^k' respectively, we find in the comparison 
of Eq. lf75|l with Eqs. and l(B7|l that we have indeed matched the terms {D^O^^) and (Dj^O^) to machine precision 
ijLiebendorfer ll2000j) . 
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3.3.5. Angular aberration from observer motion 

Motivated by this success, we try the method of characteristics also on the angular aberration corrections. In this 
case, the characteristic is described by Ct + = 0. We write the prefactor of the angular aberration correction in 
Eq. H21I) as a function of the quantities A and B defined in Eq. H70() to obtain with ( — 1 — /i^, 

0^ = iA + B/C)-^[CfiF]. 



As before, we try to eliminate the partial derivative with respec t to the angle cosine with a time derivative along the 
characteristic in 
find the relation 



characteristic in the angle part of the momentum phase space l)Liebend5rfer ll2000|) . For the quantity Ra = r^p, we 



dhiRa 
adt 



= A + B/C 



With this, we transform from the "Eulerian" variable x = fi to the "Lagrangian" variable y = C ^^^fJ-fRa- After a 
multiplication with Cm, the angular aberration correction = Ct + becomes: 



d£ 
'dt 



dlnRa d 



dt dfi 



+ C 



-1/2 dRg ^ .3/2 p d 



' d[C'/^plRa] 
dt 



d [c-^/^^l/Ra] 



dt 



The particles initially residing in the angular interval (l — 3//^) A/i = C2M2 ~ CiMi ^-^e shifted by angular aberration 
along characteristics with constant /i/ (v^i?^): 



[(1-3m^)FAm]') 



0. 



(79) 



We have to keep in mind that any correction to the particle propagation direction also affects the energy conservation 
in the frame of a distant observer. Therefore, it is desirable to construct a numerical implementation of angular 
aberration that prescribes the changes in the specific luminosity X]j Mi'^*' t° machine precision. Hence, we 
evaluate the change of the specific luminosity, dl = {\ — 3/i^) fiFAfi, along the characteristic in analogy to Eq. (|33|l . 



[(1 ~ pFAp] 



OA 
dt J 



dt 



(80) 



= (1 - 3^^) FAn 

We identify the bin size (l — 3/i|,) A/ijv = Wj' with our Gaussian quadrature weights. In the finite difference repre- 
sentation, Eq. (|79|l leads to the condition for number conservation 



Fi'.j',k'Wj' 



= 



and Eq. H80|) to a prescription for the numerical evolution of the specific luminosity 



F 



i' ,j' ,k' f-j'Wj' 



F 



— ~ iCj'Ai'^k' + Bi'^j'^k') Fi'j'^k'fJ'j'Wj'ai'dt. 

The direction of the differencing can be chosen by dj — ±1. The change in the neutrino distribution function from 
angular aberration is then: 

with 



-{Ai'^k' + Bi'j'^k'/Cj') 



'Cj' l^j' Fi' ,j> ,k'Oii'dt 



''^i',i',k' ~^i',j'-di,k'- 

Combined into one update, this leads to the following finite difference representation of the angular aberration term 
in the Boltzmann Eq. (|15|l : 



{Aii^k' + Bi'^j>-dj,k'/Cr-dj) 



Vj'-dj 



{Ai'^k' + Bifji^k'/Cj') 



fljr t'-j'-dj 

Cj' f^j' Fi\j' ,k 



Cj ' -djf-j' -djFi' ji-dj,k' 



(81) 
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We apply the aberration corrections with dj = +1 for fj. < and dj = —1 for > 0. This is not upwind differencing 
and, therefore, runs the risk of producing negative particle distribution functions. However, there are two good reasons 
to accept this shortcoming: (i) The angular aberration correction is generally small, with the exception of aberration in 
the vicinity of strong shocks with large velocity gradients, (ii) The chosen direction of dj guarantees that no particles 
are shifted off the grid. This is a prerequisite for number and energy conservation. If we calculate the contribution of 
to the energy evolution, J2j k (^i+i + Ui+i^ji) O^^Wj/E'^idEki , along the lines of Eq. 1)35(1 . this choice of dj causes 
the perfect vanishing of boundary terms in the discrete "integration by parts" with respect to the angle cosine. We 
obtain 

Ol'^ ■■ + Mj'^^'./c' - F,,j,,k'WfEldEk'. (82) 

The second and third terms in the parenthesis cancel exactly with the yet unmatched terms from 0^~^ in Eq. (|78|l . 
Thus, also the cancellations (O'^O^) and (O^O^) are guaranteed in Eq. (|37|l . However, in failed supernova simula- 
tions, the shock recedes at late time to very small radii and becomes very strong because of high infall velocities. We 
have encountered numerical fluctuations in the neutrino distribution function at the shock position at this late time 
that are due to this fixed choice of differencing in the angular advection in the aberration terms. They disturb the 
progress of the simulation with large time steps. In future long term simulations, we will use strict upwind differencing 
also in the angular aberration terms and allow a deviation from perfect energy conservation in these terms (we still 
conserve particle number). In section ^31 we demonstrate that the expected energy violations from angular aberration 
are unlikely to dominate the violations of energy conservation by the adaptive grid or a mismatch in (^CfD^D^Y 

3.3.6. Frequency shift in the gravitational potential 

The action of the gravitational potential on the propagation of the particles, Ct + De = 0, in Eq. (|19|) . is finite 
differenced in full analogy with Eq. H77() , 



De = - 



El,dEk' 



dEk'-dk p dEk' 3 

-i^k'-dk^i',j',k'-dk ~ —l^k'^i'd'M' 



i^k' — i^k'-dk i^k'+dk — i^i 



(83) 



This is a valid finite difference representation of De = —jiG^/E"^ ■ d {E^F^ jdE. We apply the corrections from Eqs. 
((77(1 and ((83(1 with dk = +1 (blueshift) if the two following conditions are met: (i) J^j k ~ Bi',ji^k') Wj'El,dEk' 

is positive, and (ii), its absolute value is larger than J2j k l^J'^i+i'^J'-^l'dEk' ■ Otherwise, we choose dk = —1 for 

a frequency redshift. With this, the chosen sign of dk in Eqs. ((77(1 and ((83(1 implements upwind differencing for the 
angle-integrated distribution function. Note, that the direction of the finite differencing in energy space must not 
depend on the individual angle bins. An angle dependence would destroy the important cancellation of redshift and 
blueshift in the limit of an isotropic distribution function. 

In Eq. 1(68(1. we have introduced the placeholder = F^, whose finite differencing we will now define. Eq. 1(41(1 
describes how the terms {Cf D^D^Oj^) should combine to the expression Airrp (1 -I- e -I- p/p) H . The latter is used in the 
equation for the total energy evolution. A matching of the energy changes in Cf, D^, O^, and to this expression is 
the last constraint of the list in Eq. 1(371) that we have not yet implemented. The term Cf is determined by a~^du/dt 
in the hydrodynamics equations. The finite differencing of the terms and Oj^ has also been determined in the 
second line of Eq. ((611) and in the first term of Eq. ((82(1 respectively. We simply resolve for and obtain 

+ 7^ V {A,,^k' - Ai',k')F,,^j,^k'lJ-j'Wj>El,dEk'. (84) 

The evaluation of the Lagrangian time derivative du/dt on the adaptive grid requires the appropriate corrections. The 
application of Eq. I(44() to the velocity leads to 

du 1 (u,da^-u,da, . ^.^j ^ roi * i \ /n.N 



dt dai \ dt 

where we have once more used the grid velocity with respect to mass, uj?' — — (a^' — a^') /dt. The advected velocity 
is evaluated by upwind differencing: 

u, if w'?' > 



Mi+i otherwise. 

3.3.7. Boundary corrections and the enforcement of Fermi statistics 

fcmax energy groups are subject to a frequency shift, but there are only fcmax— 1 enclosed group edges to shift particles 
across. In the energy phase space, upwind differencing is necessary for the stability. Therefore, the straightforward 
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application of Eqs. H77() and (|83(l leads to particles that fall aside the grid. The lost particle number is proportional to 
the particle population in the energy groups at the boundaries of the grid. At the lower boundary it is small because 
the phase space volume E^dE is small. At the upper boundary it is small because the neutrino distribution function 
is small. The particle loss at both boundaries would approach zero when the coverage of the energy grid would be 
extended. Nevertheless, we want to correct for this effect and avoid that particles and energy unaccountcdly disappear 
from the radiation field. On the one hand, we would like to support requirement (iii) in the introduction of this section 
with an accurate evolution of the radiation moments, even and especially at low energy resolution. On the other hand, 
conservation to machine precision in wellknown terms is an existential practical asset in the hunt of unknown energy 
leaks in a complex code. To resolve this issue, we simply do not apply the energy shift to the boundary group if this 
would advect particles off the grid. This exception fixes particle number conservation. The induced error in the energy 
balance of the radiation field is then restored by the application of minimal energy shifts to the population of all other 
energy groups. 

Before we describe the details of this correction, we mention another exceptional correction to the frequency shift. 
It is applied when the Fermi statistics of the radiation particles is violated by numerical effects. The collision term Cc 
includes all necessary blocking factors to avoid the overpopulation of states in the fermionic particle gas by emission or 
inscattering. If inelastic interactions occur on a very short time scale, they will efficiently restore numerical deviations 
from perfect Fermi statistics in the radiation field. However, an overpopulation of energy groups can still occur in 
the numerical implementation. Consider for example a hydrostatic situation with a coUisionless degenerate particle 
gas. The Boltzmann equation H15|l reduces in this case to Da + + De = 0. It is straightforward to demonstrate 
that pF = 1 is a solution of this equation. In the numerical implementation, however, it is unlikely that the redshift 
towards low energy groups, De, is to machine precision compensated by the space propagation term, Da^ and angular 
advection term, D^. Even if the collision term is large, overpopulated energy groups may occur as a consequence of 
operator splitting. Our previous code version l|MezzacaDna fc Bruenn lll993a() implemented the observer corrections 
in Eq. H77|) . (|81|l . and H83|) in an operator split explicit update to the implicit solution of the transport equation 
Ct + Da + Dfj^ = Cc- In this case, overpopulation is more likely to occur. The scheme described below has been 
developed for this context. For simulations with massive progenit ors in general relativ ity, however, we have expanded 
the implicit solver to include all terms of the transport equation l)Liebendorfer et al. i 2002). For cosmetic reasons we 
still strictly enforce Fermi statistics when necessary. We performed a control run with both corrections switched off. 
At any time before and including bounce, deviations in the luminosity are smaller than 2%. Differences in the rms 
energies are smaller than 0.5%, and in any hydrodynamic quantity smaller than 0.2%. When the shock propagates 
through the neutrino-degenerate core, the 2% differences in the luminosity cause transient differences of the same size 
in parts of the velocity and the density profile. Later, the differences are smaller than 1% in any quantity and at any 
time. 

As in l)MezzacaDDa fc Bruenn iflQQSaT . we recursively shift overpopulated particles into the next higher energy group. 
This procedure is number conservative, but it introduces an artificial increase of the internal energy of the radiation 
field. The details of the correction scheme for low resolution boundary effects and overpopulation are the following: 
First, we quantify the energy correction: 

Ri',j' — {fJ-'j' A-ii ^k' — Biij'^k' — A'j'G'^i) Fii ^k'ElidEk' , 

with k = fcmax, if dfc = 1 (blueshift); and fc = 1, if dfc = — 1 (redshift). The corrections stem in most cases from an 
energy redshift. Thus, we set dk = —1 and initialize the accumulation of corrections in Aiiji^k' with A = 0. Then, 
we loop over fc in ascending order and carry out the following procedure: (i) Distribute the correction Ri'j' over all 
groups £' with k + 1 < £ < fc,„ax according to 

A _ A J 

l\:t „■/ nt /_\jt At p! — — . 

Z^S = fc + l,fe,„ax ^l' ,j' ,s'-t^s'"--^s' 

(ii) Apply the correction Ai'j'_fc/+i according to the conservative finite differencing used in Eq. (|76|l 

n+ - A dEk' + l p3 p 

J^k' + l — J^k' 

and update the particle distribution function 

Fz',j',k' - Fj'j'^k' ^ "^tj^k' - ^i'.j',k'-i 

a^'dt EldEk' ■ ^ ' 

(iii) If the distribution fimction Fi/j/^k' is larger than 1/pi', correct the particle flux into this group by the amount 
of overpopulation, ji f.i = n^, ■, j,, + {1/ — Fiij/^k') E'^idEk' , and evaluate the energy change owing to this flux 
correction, 

R^',f = (IM' - i^«'j',fc') [Ek'+i - Ek') EldEk'. 

Set Fi' ji k' to the maximum value, 1/pi' and apply the energy correction Ri' j' with the next iteration of the loop over 
fc. ' ' 
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3.3.8. Collision term 

The finite difference representation of the collisio n term, Eq. H22|) , is described in detail in l|Mezzacappa fc Bruenn I 



ll993blctlMezzacaDPa fc Messer ll999HMesser ^ 20001. We represent it here with a simple placeholder, namely emissivity 
jV,fc' and opacity Xi',k' as a function of the thermodynamical state of the fluid. 



Pi' 



- X^',k'{pi'■,T*,,Yl^,)F^,J,^k' 



(87) 



The quantities with a star refer to quantities that have consistently been updated with the particle distribution 
functions before the actual hydrodynamics update. We outline the time sequence of the evaluation of equations in 
section 

3.4. Coupling between the radiation field and the fluid 

The neutrino radiation field in the supernova is tightly coupled to the high density fluid. In order to obtain a well 
defined link between the evolution of radiation quantities and the evolution of hydrodynamical quantities, we recall 
that the solution of the microscopic Boltzmann equation updates the radiation moments according to Eqs. H28|l - (|3()|l . 
We may now subtract these updates from the global evolution equations (|?^- (fT^ to find the hydrodynamical part of 
the evolution equations 
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(91) 
(92) 
(93) 
(94) 

(95) 



- - xFj E^dEdp. 

The equations are written such that the presence of the radiation field only enters in terms of energy and momen- 
tum exchange. In the limit of a decoupled radiation and matter fiow we therefore solve for ideal hydrodynamics 
and free streaming in an independent and numerically stable manner, no matter what the size of the radiation field 
is. The only remaining interactions between the radiation field and the fluid are of a gravitational nature, for ex- 
ample in the contribution of the radiation field to the common gravitational mass in Eq. (|93|l . or in the term 
47rr^ ((1 -I- e) pK + pJ) in the momentum equation (19011 . The detailed discretization of Eqs. H88|l - H95|l has been de- 
scribed in l|Liebendorfer. R.osswog feJIhiaLemflnn,. .2QQ2) . The interaction terms between the radiation field and the 
fiuid are 
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(96) 
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where the minus sign in the Ye-teim is used for electron neutrinos and the plus sign for electron antineutrinos. Emission 
and absorption of fi— and r-neutrinos do not change the electron fraction. The star superscript for the temperature 
and electron fraction variables indicates again an evaluation of the cross sections for emission and absorption with a 
consistently updated thermodynamical state as described in the next section. 

3.5. Implicit Solution 

We have now a finite difference representation for every term of the transport equation. In this section we join 
them to one implicitly determined solution. Each time step cycles through an implicitly finite differenced update of 
the neutrino distribution functions and the thermodynamic state of the matter, an implicitly finite differenced update 
of the hydrodynamics variables and the adaptive grid, and an explicitly finite differenced cosmetic correction. The 
intermediately updated values do not belong to uniquely defined physical times. We denote them by fractional values 
of the cycle counter n. In detail, the update of the primitive variables for one cycle proceeds as follows: 

1. We start with the old neutrino distribution functions F :— F{n — 1/2) for the electron flavor neutrinos. 
They are given _on a Lagrangian mass grid a. We assume that we know the old hydrodynamic variables, 
{a, r, u, m, yo, r, y"e, a} :— H{n — 1), and that the new hydrodynamic variables, {a, r, m, m, p, T, 1^6, a} '■— H{n), 
have been calculated in step (3) of the previous cycle. In one implicitly finite differenced update, we solve for 
F{n) := F as prescribed by the transport equation Ct+Da + Df^ + DE + OE + Ofj, — Cc = 0. The finite differencing 
is chosen according to Eqs. (|53)), (|57)), (77|)- EJ, (|H3), and This imphcit update additionally includes 
equations for the evolution of the temperature and electron fraction due to particle- fluid interactions. These 
quantities are strongly coupled to the radiation field. We solve for 

_^=°'-<'* 

where the source terms are given in Eq. H96|l . Note that Eq. H97() also determines the temperature because 
the internal energy and temperature in the implicit update are always consistent with the equation of state, 
e*{p,T* ,Y*) and e{p,T,Ye). We repeat step (1) for the muon flavor neutrinos. The tauon flavor neutrinos 
are included by counting the appropriate muon quantities twice. But we always distinguish neutrinos from 
antineutrinos. 

2. In this minor second update, we correct for discretization errors caused by the advection scheme in the energy 
phase space and for numerical violations of Fermi statistics. This update is not related to any physical terms 
and only relevant in low resolution simulations. We set F :— F{n), and update F(n + 1/2) :— F as described by 
Eq. jHSI)- 

3. In the update of the hydrodynamics variables, the primitives {a,f,u,in, p,T,Ye,a} := H{n) are evolved to 
H{n + 1) := {a, r, u, m, p, T, Fe, a} in a fully imphcit step according to Eqs. Note that we use the neu- 
trino distribution function F{n) rather than F{n+\/2) to evaluate the moments J, iJ, and K for the gravitational 
coupling. The former are the direct result from the implicit solution of the transport equation and represent more 
accurately the quasi-stationary radiation quantities. The hydrodynamics update is based on the transfer of inter- 
nal energy, e*"^*, electron fraction, Fg*"^*, and the neutrino stress, 5"''', as evaluated in the firs t step of this cycle. 
The detailed finite difference representation in agile has been outlined in l|Liebcndorfcr. Ross woe fc Thielemann I 

Steps 1, 2, and 3 form literally a cycle; there is no position where all data are consistently updated. After step 1, the 
neutrino distribution functions F{n) are in accurate equilibrium with the matter. After step 2, the energy conserving 
corrections from energy advection are included and the distribution functions i^(n-|- 1/2) obey Fermi statistics exactly. 
After step 3, the overall energy conservation is balanced, but the new hydrodynamic variables H{n-\-l) already live on 
a different spatial grid than the neutrino distribution function. We chose to set the integer value of n at the instance 
we dump the result into files. This is between step 1 and step 2, where all quantities live on the same grid positions 
and where all rates and neutrino distribution functions accurately reflect equilibrium conditions if the mean free paths 
are small. 

Above description is complete with respect to what we compute. It does not address the question of how we compute. 
We apply a Newton- Raphson scheme to solve the implicitly finite differenced nonlinear equations in step (1) and step 
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(3). The solution vectors for the two updates are 



. ,rad 



F" 

J"., 



m,; 



PV 

Y / 



(99) 



At the present resolution with 103 zones, 6 angular bins, and 12 energy groups, the radiation solution vector in step (1), 
y'''*'*, has the length imax x (2 x jmax x fcmax + 2) = 15038 and the hydrodynamics solution vector in step (3), y^-^'^, has 
the length (imax + 1) x 8 = 832. The radiation system is much larger than the hydrodynamics system and dominates 
the computational effort. The hydrodynamics system with the adaptive grid, however, shows a more complicated 
coupling between the variables. Common to both cases is that, in one time step, we have to determine the vector 



y 



at time = t''^ + dt which solves a nonlinear system of equations, i?(y" 



y 



n+l. 



dt) = 0. The Newton-Raphson 



scheme is based on the linearization of the system of equations around a guessed solution vector, y""*"^ , 



R(.y\y 



n --n+l 



Ay;dt) = R{y^,r^';dt) 



dRiy",r+\dt) 
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Ay + 0{{Ayf). 



The residual on the right hand side of Eq. H100|) can be neutralized by a specific choice of corrections. 
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f dR{y",y^+'-dt) 
\ dy 
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R{y'\r^';dt). 



(100) 



(101) 



~n+l 



y 



Ay is then a better guess because it fulfills 



R{y'\r^' + Ay;dt)=0 + O{{Ayy). 



This procedure is iterated until a norm of the correction vector \Ay\/y"'^'^^ becomes sufficiently small. In the solution 
of the radiation system we stop the iterations when the corrections A (pFi'ji^k'), Aei' /ei', and AY^^i' lYe,i' are smaller 
than 10~^°. If the radiation is rather stationary, (e.g. in early phases of the gravitational collapse or some time after 
bounce during neutrino heating), the solution of the transport equation requires about three iterations. Around bounce, 
the iteration count is typically 4 and may only occasionally reach 6. Convergence problems appear only at very late 
time (around a second after bounce, depending on the progenitor model) when the adaptive grid demands tiny mass 
zones (smaller than ^ 1O~'*M0) to resolve the steep density gradient developing at the surface of the protoneutron star. 
Of order half of the computation time is spent on the construction of the Jacobian 9i?(y", y"+^; (it)/9y"+^, while the 
other half is required to solve the linear system in Eq. (|101(l . We solve the linear system by direct Gauss elimination. 
The zone- wise arrangement of the solution vector in Eq. (|99|l leads to a block-tridiagonal Jacobi matrix because every 
equa tion depends o n inpu t from at most three adjacent zones. A straightforward solver for band-diagonal matrices (e.g. 
from lPress et al.1 ()1992f l or from the LAPACK linear algebra package) works perfectly if the variables and equations 
are scaled such that the maximum element in each row of the Jacobian is of order unit. A solver for band matrices 
is not only simpler than a solver for block-tridiagonal matrices, but it adapts even better to the boundary of the 
off-diagonal blocks because the latter have nonzero coefficient s only on their diagon al. Iterative solvers require careful 
preconditioning to solve the described equations satisfactorily l|D'Azevedo et al. 1 120021 . In order to calculate the exact 
Jacobi matrix, we need derivatives of the equations with respect to all primitive variables. Different methods are used 
to build the Jacobi matrix for the derivatives with respect to the particle distribution functions, the derivatives with 
respect to temperature and electron fraction, and the derivatives of the hydrodynamics equations: 

The left hand side of the Boltzmann equation H15|) is a linear operator on the particle distribution function. Thus, 
one can inquire Eqs. (|5^ . (|57l) . lf77|l . and for the occurrences of Fi'±ij'_fc', Fi>j'±i_k', and Fi>j'_k'±i- 

The corresponding prefactors, separately collected, form the contribution to the Jacobian for this part of the transport 
equation. There are no other dependencies on the left hand side of the transport equation. It requires more tedious 
efforts to gather the derivatives of the collision term. On the one hand, the interactions can in principle couple all 
angle bins and energy groups within one zone. On the other hand, blocking factors lead to nonlinear dependencies 
in some reactions. This applies as well to the derivatives of Eqs. H97|l and l|98|l with respect to particle distribution 
functions. The coefficients for the Jacobian have to be read off from extensive listin gs of sum mations over angle bins 
and energy groups, as outlined for example in Eq. (67) and (94) in (Mezzacaooa & Messer l irggQ). 

Too complicated to be handled manually are the derivatives of the collision term in Eq. (I87II . Eq. (1971) . an d II98II 
with respect to temperature and electron fraction. The dynamic table developed in ((Mezzacappa fc Bruenn ■■1993a[) 
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helps in this case. We iUustrate its mechanism with the most simple example, the derivative of the left hand side of 
Eq. 1)97(1 with respect to temperature, 

d f e*,{pi>,Ti,,Ye^,,)-ei,\ ^^^^^ 



dTi, V dt 

We discretize the three-dimensional space with dimensions p, T, Yg into cubes around the thermodynamical states of 
the fluid zones {pi' ,Tii ,Yf,,i')- The eight cube corners 

{(/?;, r„i, Fe,„) , l,m,n= 1,2} 

satisfy 

logio iP2/pi)^N-^ 
logio (T2/Ti) = 7V^i 
Y2-Y,^Ny\ 

The numbers Np, Nt, Ny define the resolution of the cube. In our applications, we use Np — 10, — 40, Ny — 50. 
If a quantity has to be evaluated as a function of pi' , T^/ , and Yg , we fetch the cube for zone i' and evaluate the 
quantity on all eight cube corners if this has not already been done in an earlier cycle. For the evaluation o f the i nternal 
energy, for example, eight calls to the interactive Lattimer-Swesty equation of state ( Lattimcr & Swcstv1 fl991t) would 
be required to store the internal energy on the cube corners. The internal energy of the zone is then found by tri-linear 
interpolation in the logarithmic internal energy, 

logio e^' = (1 - Cp) (1 - Ct) (1 - Cy) logio e (pi, Ti, Fe.l) 
+ (1 - Cp) (1 - Ct) Cy logio e (pi, Ti, i;,2) 
+ (1 - Cp) Ct (1 - Cy) logio e (pi, T2, Y,s) 
+ (1 - Cp) CtCy logio e (Pi, 1^2, ye.2) 
+ Cp (1 - Ct) (1 - Cy) logio e (P2, Ti, 
+ Cp (1 - Ct) Cy logio e {p2, T^, Y,,2) 

+ CpCT (1 - Cy) logio e (P2, 1^2, (103) 

+ CpCTCy logio e(p2,T2,ye,2) 



where 



^ logio 

logio (P2/P1) 

_ logio m^m) 

^ logio (^2/^1) 

i"e.2 - re,i ■ 



Cy 



The linearity of Eq. H103|) makes the determination of the derivatives straightforward. For example, in Eq. (|102|l we 
need 

de,> ^ e,, g logio e,- ^ .r ^ logio e^' 

OT,, T,, 9 logio T,' ^T,' dCT ^ ' 

where 9 logio ^i' /QCt can simply be read off from Eq. (|103|l . This scheme has important advantages. First, the relation 
between the values of the interpolated variable and its derivatives inside the cube is exact, both being derived from 
the same interpolation formula. This helps convergence in the multidimensional Newton-Raphson scheme. Second, 
the subroutine generating the variable on the cube corners need not be used at each time step, but only if the point 
{pii ,Tii ,Yf. ,ii) moves outside the cube, in which case a new, adjacent, cube is generated. As the complexity of the 
scheme does not depend on the complexity of the function that generates the variable on the cube corner, we apply 
the scheme to all reaction rates and scattering kernels, whose derivatives with respect to p^', T^/, and Ke.i' are then 
straightforward to determine as in the example of Eq. (|104|l . Given that only a few percent of the cubes have to be 
regenerated during one time step, these dynamic tables enhance the computational efficiency of our method by a factor 
of ~ 50. Third, we avoid discontinuous changes in the variables and any of its derivatives, which can be traumatic 
for a Newton-Raphson scheme. Note, that we occasionally allow extrapolation if the thermodynamic state leaves the 
cube during the iterations for the solution of the radiation equation in step (1). The cube boundaries are only updated 
afterwards in the hydrodynamics step (3). In the solution of the hydrodynamics equations, however, we update the 
cube boundaries for every iteration of the Newton-Raphson scheme. As the hydrodynamics equations only rely on the 
pressure and internal energy from the equation of state, this is not a computationally expensive measure, although 
it occasionally deteriorates the convergence behavior of the implicit solution. It is necessary because an extrapolated 
energy in the converged solution would not necessarily match the interpolated energy after a postponed cube update. 
This would lead to a violation of energy conservation. This problem does not arise with the internal energy in the 
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solution of the transport equation because only the energy exchange rate, e^J^*", from Eq. (|98|l is transferred to the 
hydrodynamics equations in step (3). The extrapolated energy, e*,, is discarded after step (1). 

The Jacobian of the hydrodynamics equations is much smaller. But it has a less regular structure due to the strong 
nonlinearities in the equations of general relativistic hydrodynamics. Furthermore, the adaptive grid couples changes of 
all variables to the grid motion, which itself depends on the gradients of the variables. The spatial smoothing operator 
in the adaptive grid equation couples variables from five adjacent zones. In order to eliminate the considerable risk 
of producing an erroneous Jacobian in lengthy manual derivations, the hydrodynamics code agile constructs the 
Jacobian automatically by taking numerical derivatives ijLiebendorfer. Rosswog fc Thielem ann 112 00^ . We start with 
the evaluation of the residuum vector R = R{y^'', y""^^; dt) based on an initial g uess jy""*"^ and defin e, based on machine 
precision, a small number e for the calculation of the numerical derivatives, ijPress et al. II1992D . Next, we sort the 
components of the state vector, denoted with label i, into J distinct groups gj according to the rule that none of the 
equations R may depend on more than one state vector component out of the same group. We try to choose these 
groups such that the number of distinct groups is minimized. We select a group gj and create a variation ^'""'."+1^ 

~r-ii?ar,n+l -r-in+l , _ r-in.Scl 

y[i\ ' ^ ^y[i\ ^ +ey[i\ ■ 

m='^ — M scf ' (105) 

for all components i £ gj- The residuum vector i?"'"' = R{y", ^""'"^"+1; dt) is then evaluated based on the varied guess. 
From the two residuals we can extract the components of the Jacobian 

_ ,i dR[k]{y-,r'-';dt) Rjkrr _ ^[fe] 

A[k,^\-y[^\ g^j^j„+i - ^j^j (106) 

for all i G gj. The Jacobian is complete when this procedure has been performed for all groups J. Finally, we scale 
the rows of A and the right hand side residuum vector by the maximum component in the row and solve the linear 
system 

^ A[k, ^] (AyW/y[z]"^scl^ = -i?[fc]. (107) 

i 

to get the corrections Ay for the update of the guess. Again, we use direct Gauss elimination for band-diagonal matrices. 
Before the next iteration, all zeros in the Jacobian are detected and the sparsity structure is refined. This allows further 
reduction of the number of required groups for all following time steps. The evaluation of R^"^^ — _R(?/", dt) 
for a group j is completely independent from the corresponding calculation for a different group j' ^ j. In order to 
compose a complete Jacobian, the system of equations has to be calculated once with the actual guess, and J times 
in parallel with varied guesses. The hydrodynamics solution usually converges in 3 iterations. 

The actual time step is set such that the hydrodynamics variables don't change by more than a percent per_time step. 
For the neutrino distribution functions, it has been most efficient to apply the weak condition (pF — pF) / (^pF + O.l) < 
0.1 combined with a general upper limit of the time step to 0.1 ms. Both radiation transport and hydrodynamics 
proceed with the minimum time step compatible with any of these conditions. A simulation from 100 ms before bounce 
to about 600 ms after bounce takes about 12000 time steps. The time steps are smallest around bounce, where they 
show typical values of 10"** ms. 

4. CODE VERIFICATION 

In this section we pres ent a test of the terms and finite difference representations introduced since 
l)MezzacaDDa fc Bruenn lll99 3a'). First, we investigate the diffusion limit in our finite difference representation and 
demonstrate that a small diffusive fiux is accurate in the presence of a large, nearly isotropic radiation field. We also 
check the other extreme, the evolution of the radiation moments in a free streaming situation in spherically symmet- 
ric geometry. Another test in stationary space-time investigates the implementation of gravitational terms, such as 
number luminosity conservation, gravitational frequency shift, and gravitational bending. The frequency shift and 
angular aberration of the radiation field are probed at the shock front, where we relate the discontinuity in radiation 
quantities in the comoving frame to the smooth radiation field in the view of stationary observers. Then, we investigate 
the resolution dependence of our results and perform a detailed energy and lepton number conservation analysis to 
check the overall consistency of our code. Finally, in the intended application of stellar core collapse and postbounce 
ev olution, we compare our results wi t h the independently developed multi-group flux-limited diffusion (mgfld) code 
ofH ruenn. DeNisco. fc Mezzacappal lj2001(l . The latter is based on a sophisticated but approximative treatment of 
radiative transfer in spherical symmetry and uses a different hydrodynamics code. 

4.1. Diffusion limit 

The neutrinos in the interior region of the protoneutron star are almost completely trapped. Because of the high 
electron chemical potential, they equilibrate with matter at high root mean square energies. The neutrino number 
density multiplied by the speed of light exceeds the neutrino number flux by orders of magnitude. It is therefore 
essential that the finite difference representation of the Boltzmann equation does not allow errors in the large particle 
density and pressure to swamp the markedly smaller particle flux. A derivation of the analytical diffusion equation 
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from the Boltzmann equation is the first step in the derivation of a multi-group flux-hmited diffusion scheme l|Bruenn I 
Il985j^ . We perform here exactly the same procedure, but apply all operations to the finite difference representation 
of the Boltzmann equation. For simplicity, we assume a stationary background. The particle distribution function is 
split into an isotropic and a flux component, 

We further assume that the flux component is stationary, i.e., that its time dependence in the radiation momentum 
equation is negligible. We substitute this approximation for the neutrino distribution function into the finite difference 
representation of the stationary-state Boltzmann equation, Ct + Da + D^^ + De = Cc, and obtain the radiation 
momentum equation by the application of the operator (3/2) ^ and the use of the identity (3/2) ^ ^j3,Wj> = 1: 



with 



R?, + Rn + R?, + R}, + R'p — Rc 



2-77 

- r- (Q;,/_ip,'_iV'°_i,fc- -l-a.'Pi'Vz^vfc')] 
— 37r 



- fi (1 - 2/3,^fc/) (a,/_ipi'_i7/'^,_i fc, - aepi'ipl, ,,,)] 
i?o=-2T,+i^°,,, 



J^f.,a£jk' \i^k' — -CJk'-dk i^k'+dk — -CJk' 

R\^-X^',k'{p^',T:,,YlMlM'■ (108) 

The simple expression i?" results from the specific property Ht)3|) of the angular difference coefficients. The radial 
dependence of I3i^k' in Eq. (|56|) can lead to an asymmetric treatment of ingoing and outgoing particles, producing the 
numerical term R\. Also R^^ arises as an unwanted term with no analytical correspondence. It is therefore essential 
to set the transport coefficients in the diffusive regime precisely to j3i^k' = 1/2. With the choice in Eq. H56|l we 
immediately find R\ — 0. Also R^^ vanishes because the term under the sum over j becomes antisymmetric with 
7i',fe' = A',fc' = 1/2- For the next step, we need a more concise notation for the difference and average of a zone edge 
variable, V: 

d{V)^,^V,+i-V, 

For zone center variables, we use the analogous definitions, e.g. dl^aptp^,) . = ai'pi'if}'^, j., — k'- With 

these, we may rewrite the terms i?" and i?" as 

K= \J' P^'^l^' + — — rl^,d{api:l,),^^^+rld{api,l,)^ 

^d{r% 
2d(r^ 



If the hydrodynamics equations relate the rest mass element to the volume element by Fi+idfli' — [An/3) piidlr^^^,, 
the first term in R^^ cancels with the first term in i?^. Finally we try to extract the lapse function a from the gradient 
in the remainder of R^, 

' {r'd{api.t,))^,^-l-{Anr^{pr,,)da)^^, + -^ 



and find for the diffusion limit in the finite difference representation 

-Ka + + -Kb - R 

with 

K" = ^{4nr'ad{p^jl,))^, 
ai'dai' * 
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po _ '^^li dEk'-dk 3 ,0 dEk' 3 

^-E^ Ei2 JIT I 1? 7? ^k'-dkri',k'-dk T^-'^k'Vi' ,k' 

EydEk' \Ek'~Ek'-dk Ek'+dk~Ek' 

Rl^-X^'M{P^',T*,,Y:,,)^l^l,,. (110) 
If we compare this to the analytic radiation momentum equation in the diffusion hmit, 

where we have set J — 3K = 0, we can identify the physical meaning of all remaining terms. The diffusive flux 
appearing in is determined by the opacity x and operators acting on the particle density. The term in R'^ describes 
the flux generated by a gradient in the particle density. It is the only term that survives in the Newtonian limit. In 
full general relativity, we have an additional gravitational force acting on the radiation particles. It contributes to the 
particle flux and is described in the term R'^. Its analytical analogue is R'^^ = 3 (F/a) {da/dr) ijp . The factor of three 
stems from our choice 3^^*^ = "iK = J. 

The term i?^ corrects for the fact that the gradient in R'^ has been taken at constant local particle energy instead 
of constant particle energy at infinity. The terms R'^ and R\ emerge in a perfect way from our finite difference 
representation because of the special measures taken in the finite difference representation of the propagation terms 
Da and I?^ in (Mezzacappa & Bruenn 1993a). However, the requirement of global energy conservation ties the finite 
difference representation of this basic choice to many other terms, as outlined in table (|37(l . These constraints impose 
a rather obscure finite differencing in the diffusive limit for the general relativistic corrections i?^ and i?^. However, 
each term is recognizable as a finite difference representation of a physical term. If we integrate the corrections over 
energy, the contributions from i?^ are guaranteed to vanish. We therefore focus on the remaining correction i?^. 
In order to ensure a vanishing flux in equilibrium, it remains to be shown that the results from our finite difference 
representation do not significantly deviate from a natural finite difference representation of Eq. Hlll|l , which becomes 
after an integration over energy, 

^5^^("Vi^")^-xH» (112, 

To this purpose, we collect from Eq. H110() the finite difference representations 

-^dr.i?;° = ^— ^ ( Wad 
Ti+i 47r [rf^^ + r^+iTj + rf) at' 

-dr,,i?'0 = Unr^ (p^P) da) ., 



Ti+i ^ 47r [rf^^ + r^+ir^ + rf 

+ (r,;+i-rOGf+i2^ (113) 

and compare their sum with the natural finite difference representation of Eq. pi2|) . 

^d{a'p^l,)^. (114) 

The comparison is performed in a time slice in the evolution of the 40 M© progenitor of ijWooslev fc Weaver1ll995|) at 
400 ms after bounce. We force the hydrodynamics to be static and evaluate the finite difference expressions (|113|l and 
H114|l in the high density region, where the diffusion limit is valid and gravitational terms sizeable. We find in Figure 
l|l()(l that the terms stemming from R'^ (dashed line) and i?^ (dash-dotted line) compose a sum (solid line) that agrees 
well with the more accurate finite difference representation given in Eq. (|114|l . Thus, we conclude that the opacities 
safely determine the diffusive flux in our finite difference representation of the Boltzmann equation. 

4.2. Redshift, gravitational bending, and the evolution of angular moments 

In this subsection we test free streaming in spherically symmetric geometry. This is the opposite limit to the diffusion 
investigated in the previous subsection. For the time being, we stay with our stationary time slice at 400 ms after 
bounce, but focus on radii larger than 40 km. Outside this radius, we switch all interactions between the radiation field 
and matter off. This facilitates the comparison of the neutrino density, neutrino number flux, and neutrino luminosity 
with the analytical behavior. Unlike in the Newtonian limit, they are subject to time lapse, gravitational frequency 
shift, and gravitational bending. 

For a stationary neutrino flux in the free streaming limit, Eq. H24|l expresses particle number conservation 

4- Unr^apH^] ^ 0. (115) 
oa ■' 

Figure Hll(l shows the radial dependence of the locally observed neutrino number luminosity Aur^pH'^ (dashed line) 
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Fig. 10. — The sum of the terms R'^ and i?,^ in Eq. IllOi determines the diffusive flux in the discrete Boltzmann equation. We check 
these terms in a stationary-state situation at 400 ms after bounce in the 40 Mq model. Plotted are the contributions from R'^ (dashed 
line) and R'^ (dash-dotted line) as determined in Eq. Illiji from our finite differencing of the Boltzmann equation. Their sum (solid line) 
is compared with a straightforward finite difference representation of the same quantity (circles) according to Eq. 11141 (The circle at 5 
km radius is the zone edge of our innermost zone, the other quantities are evaluated on zone centers). The agreement demonstrates that 
we obtain an accurate diffusion flux in the numerical solution of the Boltzmann equation. 



and the conserved number luminosity Anr^apH^ (solid line). We find that the latter is constant to machine precision 
as required by Eq. H115|) . This is the direct result of our finite differencing of the term Da, because we have placed the 
lapse function a inside the space derivative. The lapse function in the bracket converts the locally advected particle 
number per proper time to the particle advection per global coordinate time. There is a similar equation for the 
conservation of the energy luminosity. Eq. (|30(l for the particle energy carries a second lapse function inside the spatial 
derivative from the gravitational frequency shift, 

d 

— Unr^a^pH] = 0. (116) 
oa 

Figure H12|l shows the radial dependence of the locally observed luminosity Anr^pH (dashed line), the luminosity 
with the time lapse correction only, Anr^apH (dash-dotted line), and the conserved luminosity, Anr^a^ pH , with the 
gravitational redshift included. We see again, that the latter is fairly constant in the interaction free region. Unlike 
the number luminosity, the conservation of energy luminosity is not automatically fulfilled by the finite difference 
representation, and small deviations can be observed. As the local mean energies of the moving particles are determined 
by the ratio of energy and number flux, we may conclude from the fulfillment of Eqs. I|115|l and Ijlltill that, in this 
regime, the gravitational energy shift of the particles is accurately implemented. 

In order to fully constrain the redshift, gravitational bending, and evolution of the angular moments, we need to 
test the radial dependence of a third quantity. Most appropriate is the particle number density. However, some special 
care is necessary to clearly separate the following three effects on the angular neutrino distribution: (i) gravitational 
bending, (it) the changing opening angle as a function of the geometric distance from the source, (iii) the numerical 
bending caused by limited angular resolution in our finite difference method. We find an analytical expression for the 
radial dependence of the particle density by applying the operator J dEdpE^ p,~^ to the interaction- free stationary 
state limit of the Boltzmann equation (|15() , 

— {ATTr^apJ^)+T[ — 

aoa \r a or 

,2 



[j F{n = 0)E^dE - J [F - F{p = 0)] E'^dE ^ dp^ = 0. (117) 
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Fig. 11. — Shown are the radial variations in the locally observed neutrino number luminosity Airr^pH^ (dashed line) and the conserved 
number luminosity 4nr'^apH^ (solid line) in a stationary-state situation at 400 ms after bounce in the 40 model. Eq. 11151 is fuUfiUed 
by construction. 



As before, we may compare this to our finite difference representation. The apphcation of the operator 
^dEk'Wj'E^ifiJ,^ to the finite difference representation of the Boltzmann equation leads to 

(47rr-+ia,+ip.j+i J^i - Airrfa^p^J^) = 
-ai>Ti+iy^ Fi> ^j^k'dtti' El,dEk' {fif -fij'-i). (118) 

Here, atpiJf is interpolated according to Eq. H58() and Fi'j,k' according to Eq. H64() . Figure (|13|l compares the 
left hand side of Eq. H118|l (circles) to the right hand side (solid line), evaluated in our 400 ms post bounce time 
slice. The match is to machine precision; it demonstrates that we made no errors in the derivation of Eq. (|118|l or 
its implementation. We split this quantity into the contribution from geometric changes according to the opening 
angle of the source, proportional to T/r, (dashed line) and the gravitational bending, proportional to (F/a) (da/dr), 
(dash-dotted line). Because the latter is much smaller, we exclude the gravitational bending contribution to T as a 
primary source of uncertainties. We can then focus on the comparison of the analytical integral in Eq. (|117|l with its 
finite difference representation in Eq. H118() . We start with a more accurate numerical evaluation of the integral in Eq. 
H117(l in order to compare with Eq. H118|l in our time slice. The application of a Gaussian quadrature in Eq. H117|l 
converts the integral to the sum 

- a,.T,+i (^Aj2F{0),,^k'EldEk' - ^ [F,,.,.,^, - F{0),,^k']da,,EldEk'^-^^w,^ (119) 

The result of Eq. H119|) depends quite sensitively on the choice of the distribution function in the tangential direction, 
F (0) j, f., . We determine it here by the maximum entropy model for Maxwell-Boltzmann statistics in angle space, 

= Aexp(B/ij/), (120) 

(see e.g. l|Smit. van den Horn, fc BludmajTI 12000(1 for an overview on maximum entropy closures). The zeroth and 
second moments of Eq. p2Uf) define the flux factor, h — H/ J, as a function of the parameters A and B,^ 

2A 

J= — sinh(B) 
B 

* Note that our definition of moments in Eq. I23i is twice as large as in standard references because this leads to more consistency 
between our code and its description. 



42 



X 10 



53 




80 90100 



300 400 



Radius [km] 



Fig. 12. — Shown are the energy luminosities ^ttt^ pH (dashed Une), Airr^apH (dash-dotted Une), and iix r^ pH (soHd Une) as functions 
of radius in a stationary-state situation at 400 ms after bounce in the 40 model. Deviations from Eq. 11161 do not exceed 1.5% in the 
third expression for the luminosity, which is expected to be conserved. The other two expressions illustrate that the gravitational effects 
contribute substantially in this test. 



/i = coth(B) - (121) 

We derive the free parameters A and B by the inversion of Eq. (I121|l from the zeroth and first angular moments of 
the numerically obtained neutrino distribution function and define F{Q)ii^k' = Aii^k'- The result is shown with cross 
markers in Fig. H13I) . We observe comparable values between the two finite difference representations 1118() and (|119|l 
at smaller radius, where the fiux factor is smaller (~ 0.75). Large differences are found at larger radii, where the flux 
factor would be expected to be close to one. I.e., the primary source of inaccuracy in the particle density stems from 
the fact that the sum over angles on the right hand side of the finite difference Eq. IjllSI) is a poor representation of 
the angle integral in the analytic Eq. 1)1171) . 

In order to study this discrepancy in more detail, we construct a series of analytical particle distribution functions 
according to the maximum entropy model in Eq. (|120|1 with B — {0, 1, 2, 5, 10}. This corresponds to flux factors of 
h = {0, 0.31, 0.54, 0.8, 0.9}. Then, we evaluate the integral on the right hand side of Eq. I|118|) in different ways. First, 
as given by the finite differencing in Eq. I|118|l itself, then, with the more natural Gauss quadrature in Eq. (|119|l . and 
finally by an analytic integration. The analytic integration yields 



1 I d_ 

B^ ( B^ B^ B'^ 



-2-^^+^ 1^+3:3! + 5:5! + 7:7! (^22) 

The result for different flux factors and angular resolutions is shown in Table lO- We find quite poor convergence for 
the evaluation of the integral as it emerges from the finite differencing in the Boltzmann equation. The comparison 
with the center block, where no upwind differencing was used for the advection terms, demonstrates that the advective 
terms play a major role in the accuracy of the representation of this integral. A forward peaked radiation field shows 
a steep gradient in the angular distribution. The choice of upwind differencing always underestimates the advective 
flux at the edge of the angular zone. This is necessary for the stability of the scheme for arbitrary large time steps, but 
can lead to a substantial underestimation of the integral in Eq. H117|) . The most populated forward bin, for example, 
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Fig. 13. — This figure shows radial increments of the particle number density (actually plotted is the difference in the particle density 
across a zone, multiplied by Anr^). The particle number luminosity is constant in the free streaming limit and the changes in the particle 
density depends on the angular distribution of the radiation field. The dashed line shows the changes induced by the increasing distance 
to the neutrino source in Newtonian geome try. T he dash-dotted line shows the contribution from gravitational bending. The sum (solid 
line) represents the right hand side ol^ Eq. 11181 . It exactly matches the left hand side (circles). Thi s indicates a correct derivation and 
implementation. However, the crosses mark a more accurate solution to the rig ht ha nd side of Eq. 11171 than the finite differencing imposed 
by the discrete Boltzmann equation. It has been evaluated according to Eq. 11191 . The discrepancy is moderate at flux factors h < 0.75 
and large far from the source where h should converge to 1. 



never contributes to the integral. The distribution function asymptoticaUy approaches a highly populated forward bin 
while all other angular bins are emptied. The maximum achievable flux factor is smaller than one. Although higher 
angular resolution improves the situation, it does not eliminate this systematic effect. The choice of a higher order 
advection scheme in angular space might be pro mising. More rigorous, however , appears the implementation of an 
adaptive angular grid in the transparent regime (Yam ada. Janka. fc Suzuki Ill999t ). It would be designed to minimize 
the advective flow in angle space such that the undesired effects are eliminated at the root. This difficulty with angular 
advection does not occur in methods that solve a model Boltzmann equatio n along particle characteristics to close 
the s ystem of radiation moments equations with a variable Eddington factor l|Burrows et al. l[2000t iRampp fc Janka I 
120021) . An independent and much smaller source of inaccuracy stems from the fact that the integral in Eq. I|118|l is 
not written as a sum over a function evaluated at Gaussian quadrature points and multiplied with Gaussian weights. 
Hence, it does not take profit from the fast convergence of Gaussian quadrature. This can be seen if the result is 
compared with the evaluation based on the Gaussian quadrature as shown in the third block in Table J^J- 

In our supernova application, however, stationary-state convergence test s have shown that 6 angle Gaussian quadra- 
ture produces physically reasonable results |Mcsscr ct al. 1998; Mcsscrl bOOOj) . We confirm this in section 1131 in a 
comparison with a dynamical run discretized with 12 angular bins. In the regions of free streaming, where numerical 
angular diffusion becomes apparent, the neutrino field is already decoupled from the matter and does not influence 
the dynamical evolution anymore. 

4.3. Observer corrections 

For the check of the observer corrections, we go back to the original solution at 400 ms after bounce. Both, the 
neutrino luminosities and rms energies show a discontinuity across the shock front because of the Doppler effect and 
angular aberration. In the laboratory frame of a distant observer, however, the radiation field should not be affected by 
the shock. We may therefore test the implementation of the observer corrections by transforming the electron neutrino 
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Table 3 

Integral on the right hand side of Eq. (I117II . 





B 


flux factor " 


J m a,x — 6 


J ma,x — 1 2 


jmax — 48 


analytical 


Eq. JTTHJ 








-2.00 


-2.00 


-2.00 


-2.00 




1 


0.31 


-1.10 


-1.23 


-1.36 


-1.41 




2 


0.54 


-0.21 


-0.22 


-0.28 


-0.32 




5 


0.80 


0.28 


0.47 


0.68 


0.77 




10 


0.90 


0.06 


0.14 


0.26 


0.31 


Eq. liisj, 








-2.00 


-2.00 


-2.00 


-2.00 


but no 


1 


0.31 


-1.48 


-1.43 


-1.41 


-1.41 


upwind 


2 


0.54 


-0.47 


-0.36 


-0.32 


-0.32 


diflf. 


5 


0.80 


0.66 


0.74 


0.77 


0.77 




10 


0.90 


0.27 


0.30 


0.31 


0.31 


Eq. ITT^ 








-2.00 


-2.00 


-2.00 


-2.00 




1 


0.31 


-1.41 


-1.41 


-1.41 


-1.41 




2 


0.54 


-0.32 


-0.32 


-0.32 


-0.32 




5 


0.80 


0.77 


0.77 


0.77 


0.77 




10 


0.90 


0.32 


0.31 


0.31 


0.31 



"This table investigates the accuracy of the integral on the right hand side of Eq. 11171 in a parameter study based on analytic distribution 
functions with different flux factors. The result of the numerical integrations are give n for discretizations with jmax = 6, jmax = 12, and 
imax = 48 angular bins. Th e uppermost block shows the evaluation according to Eq. 11181 . It emerges from the chosen finite differencing 
of the term in Eq. I62i and converges only slowly to the exact solution. The center block shows the same evaluation without upwind 
differencing for the advected F (i.e. with 7j/ y = 0.5). The convergence is much better in this case. The lowermost block shows the 
evaluation based on the more natural Gauss quadrature 11191 . which converges very rapidly to the exact result. 



luminosity and rms energy from our comoving frame to the rest frame. In l)Liebend5rfer. Mezzacappa. fc Thielemann I 
1200 ID we determined for the particle propagation angle, ^.^ , and particle energy, , in Schwarzschild coordinates the 
relationship 

(1-^^^) {T^)-' = {l-^,^){T + u^,)-' 

r^E'^ = {T + uti)E, (123) 

where F"^ = (1 — 2m/r)^^^. Assuming free outwards streaming with transport coefficients /3 = 1 in the interesting 
region around the shock, we evaluate for the luminosity, , and rms energy, {E^"^, the finite difference expression 



Lf = 4^r2 J2 P^'-lF^'-l,J,kEf^J, ,k' El dEk' 4 ,r cwy 



(^E(^5sfcO'^«',j'.fc'^''^^fe'^.'j \^F,,y,k'EldEu'Wy]^ . (124) 

The result is shown in Figure (|14|l . The solid line is the luminosity in the rest frame according to Eq. H124() . The 
dotted line shows the luminosity in the comoving frame for comparison. We find that the luminosity profile shows a 
numerical local distortion on top of the shock front. However, if we compare the left hand side values with the right 
hand side values, we find the expected disappearance of the discontinuity, indicating that the observer corrections in 
the comoving frame are well represented. The same behavior can be seen in the rms energies. 

4.4. Radiation pulse propagation 

We have now performed many checks of the particle transport in stationary-state situations. In this subsection, we 
check the dynamics of the radiation field itself. The neutrino transport in the supernova is most dynamic when the 
shock crosses the neutrinosphere. As soon as the shock has propagated to densities where the opacities are low enough 
that neutrinos can escape from the hot material behind the shock, immediate electron capture will occur to refill the 
emptied phase space with new electron neutrinos. A neutrino burst carries an energy of order 10^^ erg throughout 
the star towards the distant observer. Although the neutrino burst is not completely interaction free, it is an excellent 
example for a time-dependent radiation field where pulse propagation can be studied. We do this by plotting the 
luminosity profiles at selected times after the launch of the neutrino burst. We shift the radial coordinate, r, in each 
time slice by the amount a free massless neutrino would have propagated during the time t since bounce, r' = r — ct. 
With an ideal numerical solution to free propagation (and assuming a point source), we would expect congruent pulse 
shapes from each time slice at the same positions r' . Figure l|15|) shows the neutrino burst with the described radius 
adjustment in the evolution of the 40 M© progenitor. We first observe, that the position of the pulse is fairly stationary. 
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Fig. 14. — Shown are the electron neutrino luminosity (solid line) and rms energy (dashed) line after the transformation from the 
comoving frame to the rest frame of a distant observer in the evolution of the 40 Mq model at 400 ms after bounce. Aside of a local bump, 
the two profiles show the expected continuous transition at the shock front. Also shown is the luminosity in the comoving frame (dotted 
line) to illustrate the effect of the transformation. 



Therefore, the pulse indeed travels with speed of light through the star. However, we also observe that the shape of the 
pulse broadens considerably with ongoing time. From the check of the neutrino number luminosity conservation in Fig. 
1)11(1 . we already know that the number and energy of the emitted neutrinos is conserved during their propagation to 
larger radii. The broadening is li kely to occur by artificial diffu sion in our numerical finite differencing scheme. Indeed, 
from considerations in iLiebcndo rfer. Rosswog &: Thielemann |[20 02:') we remember that first order upwind differencing 
introduces a numerical diffusivity proportional to the advection speed and zone width, 

A'=cdr,-. (125) 

We quantify the influence of this effect to the evolution of our neutrino pulse by approximating the diffusion between the 
time steps with AL ~ D {d^L/dr^^ At. We evaluate the diffusivity p25|l and the second derivative of the luminosity 
at the pulse peak in the first time slice. The time span At to the next time slice gives us an estimate AL for the change 
of the peak luminosity. We subtract AL from the pulse height and draw a horizontal line on the right hand side in 
Fig. H15() at this estimated peak level for the next time slice. The close to perfect agreement with the numerically 
evolved pulse height leads to the conclusion that the diffusivity (|125|l introduced with the low order advection scheme 
fully accounts for the observed pulse broadening. The resolution study in the next section illustrates that this is the 
dominant effect. 

4.5. Resolution 

In spite of our ambition to produce accurate results already at low resolution, the reader has been left wondering 
about the actual resolution dependence of our data. In this subsection, we demonstrate reasonable convergence of 
our results by separately doubling the resolution in each phase space dimension in otherwise identical simulations. 
Historically, our standard resolution with 103 adaptive zones, 6 angular bins, and 12 energy groups formed as the 
minimum resolution where we felt safe about our physical conclusions. In Fig. H16|l we compare the results with a 
simulation that uses 206 adaptive zones, one that uses 12 angular bins, and one that uses 24 energy groups. We focus 
on the quantities and instances where we find the largest deviations. We start with two limitations that are well-known. 
Graph (a) shows the mean flux factor at 100 ms after bounce in the neutrino heating phase. It can clearly be seen, that 
the run with 12 angular bins determin es the flux factor with more accuracy in the free streaming limit (Messer et al. I 
119981 lYamada. Janka. fc Suzuki Ill999ri . The angular resolution determines how close to the radial direction the most 



46 



35 



30 



IF 25 



0) 

(M 

Ln 
O 



CO 

o 



E 
3 



20 



15 



10 





-10000 



6.0ms 
9.5ms 
13.5ms 
18.4ms 
24.2ms 
34.5ms 




-5000 



5000 



Radius [km] 



Fig. 15. — Shown is the radial profile of the electron neutrino luminosity at different times when the neutrino burst propagates through 
the star (40 Mq model). The radial coordinate of each time slice has been shifted to the left by the distance a signal with light speed would 
have traveled since bounce. This makes the pulse profile roughly stationary in position. The levels on the right hand side in the graph give 
estimates of the peak height for the next time slice. The evaluation is based on the estimated numerical diffusion from first order upwind 
differencing in the free streaming regime in the transport equation. The expected peak heights accurately explain the visible decay in the 
luminosity pulse. 



forward peaked angular bins are. The closer they are the larger is the asymptotic limit of the flux factor. In graph 
(b), we show another previously documented effect |MezzacaDDa & Bruenn 19933). ^^e last half of a millisecond 
before bounce, insufficient resolution of the energy phase space leads to a poor representation of the Fermi energy 
in the degenerate neutrino gas when the trapped neutrinos are compressed to high densities. The consequence are 
rapid and transient local displacements of neutrinos. At poor energy resolution, the Fermi energies rather increase in 
steps than continuously, and the neutrino fluxes try to balance the numerical variations in the Fermi energies between 
adjacent zones. Therefore, the effect is best seen in the electron neutrino luminosities as shown in graph (b). While 
the conserved lepton number cannot show numerical variations, the electron fraction and neutrino abundances also 
reflect the transient steps in the Fermi energies. Variations in these quantities are of order 5%. However, because the 
neutrinos are trapped and because our scheme is conservative, this transient wiggles do not lead to any differences 
that survive bounce. 

This is shown in the two graphs (c) and (d) at bounce where the entropy and luminosity profile, respectively, have 
again converged with respect to energy resolution. At this stage, we detect an influence of the spatial resolution. 
During collapse, there is no special region with a high concentration of grid points. At bounce, however, the grid 
points speed inwards to resolve the newborn shock wave. We detect an effect of this rapid grid displacement in the 
entropy and electron neutrino luminosity profiles to the extent shown in graphs (c) and (d). Differences in the entropies 
can also be seen at the interface between the silicon layer and nuclear statistical equilibrium at a radius ^ 1000 km. 
The difference is of no concern because it stems from the granular triggering of zone conversions from silicon to NSE. 
Each run determines autonomously when output files are dumped. For a given time, we then compare the output with 
the closest available output files in other runs. Hence, it can happen that the conversion of a zone took already place 
in one run, but not in another one. The location of the transition to NSE shows an uncertainty of at least one zone 
width. The first 10 ms after bounce are probably the most dynamical time in the simulations. At this time we find 
again the most prominent resolution dependencies in the entropy and luminosity profiles. It is still a dependency on 
the space resolution alone. The entropy in the high energy resolution run carries the slight enhancement from before, 
overlapped by a slightly narrower and deeper cooling at the launch of the neutrino burst. The higher luminosity peak 
in graph (f ) at higher spatial resolution comes not unexpectedly, as we have seen in subsection l4.4l that the decay of the 
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Fig. 16. — Resolution dependence of the results. The standard run (circles) was calculated with 103 adaptive zones, 6 angular bins, and 
12 energy groups. We compare with runs with 206 zones (solid line), 12 angular bins (dash-dotted line), and 24 energy groups (dashed 
line). 
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Fig. 17. — The diflferent components of the energy budget are shown. The internal energy of the fluid is represented with a thin solid 
line in the upper half of the figure. We set the free constant in the internal energy such that we start with zero total energy (thick solid 
line). Most of the internal energy is balanced by the gravitational energy (thin solid line in the lower half of the figure). We show also the 
kinetic energy of the fiuid (thin dashed line) and the surface work exerted at the border of the computational domain (dotted line). The 
dash-dotted line represents the energy of the neutrinos in the computational domain and the dashed thick line represents the accumulated 
energy of the escaped neutrinos. The total energy is nicely conserved during core collapse and exhibits a systematic increase of order 10% 
of an explosion energy during the crucial phase around 100 ms after bounce. 



outwards propagating neutrino burst depends on the zone width. We can also demonstrate by this resolution study, 
that the described pulse spreading does only marginally depend on the angular resolution. Note that the apparent 
difference in the run with higher energy resolution is not a real difference, it stems from an insufficient time match 
between the output files and shows the rapidly decaying luminosity profile at a slightly earlier time. 

Finally, we compare the runs with different resolutions during the important neutrino heating phase, e.g. at 100 ms 
after bounce. Graphs (g) and (h) show again the entropy and luminosity profiles. We are happy to report that the 
simulations are converged at this stage. With the exception of the slightly higher entropy in the outer layers of the 
high space resolution run, none of the previously discussed differences has survived to this time. Moreover, the quality 
of agreement presented in these graphs is similar to what we find in the velocity, logarithmic density, electron fraction, 
and rms neutrino energy profiles at any time during the simulations. We conclude that convergence issues are far 
from affecting our physical conclusions from the simulations. Only if one asks for high precision numbers in specific 
quantities, we may, with descending importance, recommend an increase of space-, energy-, and angle-resolution. 

4.6. Energy conservation 

After having put much effort into the consistent finite differencing of the transport equation in favor of an accurate 
evolution of expectation values, we investigate in this section the energy conservation properties in the most challenging 
run started from the massive 40 Mq progenitor star. Our scheme preserves lepton numbers to machine precision by 
construction. In Fig. H17|l we show how the total energy in the simulation divides up into different energy forms during 
the simulation. The gravitational energy (thin solid line in the lower part of the figure) and the internal energy of the 
fluid (thin solid line in the upper part of the figure) form the largest contributions. For this figure, we have set the 
undetermined constant in the internal energy such that the total energy vanishes at the start of the simulation. The 
energy stored in neutrinos (dash-dotted line) grows with respect to the internal fluid energy until bounce. Afterwards, 
it evolves almost proportionally to the internal fluid energy. The kinetic energy (thin dashed line) also grows during 
collapse. It peaks at bounce and settles on a lower level afterwards, slightly decaying during the stationary postbounce 
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Fig. 18. — For the six chains of terms that are supposed to cancel in Eq. i37i we show the size of all positive contributions to 
measure the importance of accurate cancellations. The implementation of the 0(i>/c) cancellations (D^O^) and (f^O^) are most 
important (thick solid line and thick dashed line, respectively). Of similar importance is a matching in the general relativistic term 
{Cf D^D^Oji^) — 47rr(l + e + p/p)H (thin solid line). About one order of magnitude less important is the matching (C^D^^D^) (thin 
dashed line). Of lowest importance is the matching among observer corrections themselves in (0|,0^) and (0|?0j^'') (dash-dotted and 
dotted lines, respectively). The figure illustrates the steep increase in the challenge to conserve energy after bounce. 



phase because of the decreasing density of infaUing materiaL The dotted hne at the center of the figure is the work 
exerted on the surface of the computational domain. The thick dashed hne is the accumulated energy emitted by 
neutrinos. Its steep increase around bounce is delayed with respect to the other energy contributions because of the 
delay in the neutrino burst and the propagation time to the surface of the computational domain at 10* km radius. 
The thick solid line represents the evolution of the total energy in our simulation. It is very accurately conserved 
during the core collapse phase. It shows a perturbation of order 5 x 10*^ erg in the most dynamic phase around bounce 
when the grid points rush to the center to resolve the shock front. It systematically increases afterwards and reaches 
the order of an explosion energy at the end of the simulation when the neutron star collapses to a black hole. We 
further investigate this energy violation in the next two figures. 

We would obtain energy conservation to machine precision if we could enforce perfect cancellation in the six chains 
discussed after Eq. ^3, i-C (0^0^^). [D^^Of), (OlO^), {Of Of), {CjOfDl), and {CtDlD],Ol) ~ ^i^r\\ + e + 
pI p)H. The hydrodynamics scheme is designed to absorb the energy exchange by the collision integral to machine 
precision and conserve energy perfectly, even on the adaptive grid. By selecting only the positive contributions in 
the canceling terms we obtain a measure of the importance for an accurate cancellation. After an integration over 
the computational domain, Fig. H18|) illustrates this measure for the six above-mentioned expressions with a thick 
dashed line, a thick solid line, a dash-dotted line, a dotted line, a thin dashed line, and a thin solid line, respectively. 
After bounce, the maximum individual contribution to the energy conservation equation in Eq. I|27|) reaches a typical 
level of 5 X 10^^ erg/s. Fig. (|18|l makes immediately evident that maintaining accurate energy conservation is more 
challenging after bounce than before bounce. Moreover, we note that the 0{v/c) terms (D^O]j) and {DjfO^) are 
much larger than the higher order terms {0^0^} and (OfOf). 

Nevertheless, these terms are perfectly matched in our implementation. Violation of energy conservation therefore 
stems from the terms {CfD^D"^) (the matching is tuned for an isotropic neutrino distribution) and {Cf D'^D\,0^^) — 
47rr(l + e + p/p)H (the matching is based on an energy fiux averaging) as discussed in the context of Eqs. (|68|l and 
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H84|) respectively. In section IS?^ we have discussed energy conservation violations by the adaptive grid corrections when 
they are applied to the radiation quantities. Where do they enter the conservation check? If we evaluate Eq. (|27|l on 
the adaptive grid according to the recipe in Eq. (|44|l . we note that the integration of the energy over the whole star 
reduces the adaptive grid corrections to surface terms. These surface terms vanish because the grid velocity is zero at 
the center and the surface of the star. If we compare with the more detailed Eq. (|37|l we find that this time derivative 
corresponds exactly to the term C} + Cf. The energy violations by the adaptive grid show only up in the terms 
and Cf. If these terms are evaluated with all the grid corrections in the time evolution of F, u, J, and H, they will 
numerically differ from the terms and Cf we have used in Eqs. (|69|l and (|84|l for the approximate matching. We 
check the influence of the adaptive grid corrections by a comparison with a run using a pure Lagrangian grid, where 
no adaptive grid corrections can compromise energy conservation. However, in order to get enough resolution in the 
run with the fixed grid, we had to run with 400 spatial zones instead of the 103 we used with the adaptive grid. 
This run is extremely slow. On the one hand, the solution vector is four times larger. But much more important is 
that every single zone has to change its value from the preshock conditions to the postshock conditions in the allowed 
1%-change steps. This requires an almost 10 times smaller time step than with the adaptive grid, where a zone can 
follow the shock. In the latter case, the conditions in the zone changes on a much longer time scale determined by 
the drift between the zone speed and the shock propagation. In order to let a run reach the interesting phase around 
100 ms after bounce in reasonable time, we had to reduce the angular resolution to only two angular bins. Both 
measures, the increase of spatial and time resolution and the decrease of angular resolution can in principle affect 
energy conservation. Nevertheless, we hope to get the correct impression of the influence of the adaptive grid on the 
energy conservation. The two terms Cf and Cf of the Lagrangian run are also shown in Fig. H19|l (thin lines). We find 
that the order of magnitude of energy violation in the cancellation (CfD^'^D^) does not significantly change on the 
fixed grid. However, the cancellation in the term (CfD'^D^Oj^) — 47rr(l + e +p/p)H is greatly reduced on the fixed 
grid. This suggests that the adaptive grid corrections of the radiation quantities are the dominant remaining sources 
of energy violation, about five times larger than the mismatch in (CfD^'^D^). 

The absolute necessity of accurate energy conservation may be discussed. But certainly, it provided an invaluable 
check for the congruence between the programmers intention and the actual implementation of the many intricate 
finite difference expressions in our code. 

4.7. Comparison with Multi-Group Flux- Limited Diffusion 

Following our tradition ijMezzacanpa fc Br uenn 1993al), we conclude this paper with a comparison with simu- 
lations using the Multi-Group Flux-Limited Diffusion (mgfld) approximation. We note that, excepting the pi- 
oneering Boltzmann solver of l| Wilson 1 119711^ . codes with MGFLD c urrently provide t he only alternative numeri- 
cal data for the evolution of a supernova in full general relativity l)Mvra et al I Il987t iSchin der fc Bludman ' '198S 
ruenn, DcNisco. & Mczzacaooa 2001). We compare the evolution of a 13 progenitor (Nomoto & Hashimot(^ 
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1988l|) . This is a stellar model with a small iron core which — in the hope of seeing an explosion in numeri- 
cal simulations — has been used throughout the supernova literature. Questions concerning the accuracy of the 
MGFLD approximation in the dynamic semi-transparent region between the neutrino sphere and the heating re- 
gion ha ve been posed and answered in several studies with different flux limiters (in the supernova context e.g. 
( j[gnk^_^92; Mcsser ct al. 1998; Yamada. Janka. & Suzuki 1999)). For general relativistic mgfld simulations 



Bruenn. DeNisco. fc Mezzacanoa Il2001f) Bruenn developed a new flux limitcr that consists of two parts ( B ruenn 1 



l200a : The first part is a specific implementation of the usual scheme for interpolating between the optically thick 
diffusion regime and the optically thin free streaming regime, namely 

•^intrp {E) ^ 1^1 + -A, {E) ^o(^) j ' 

where At {E) is the total energy-dependent transport mean free path. This ensures that 

3 or 

in the diffusion limit, and that -ip^ (E) = ijp (E) in the free streaming limit. This part alone suffers from the generic 
problem of too rapid a transition to the free streaming limit when matter goes from optically thick to optically thin 
abruptly. To avoid this problem, a second piece of the flux limiter is constructed. It basically prevents the neutrino 
angular distribution from becoming more forward peaked than the geometrical limit. At radius r, the second part 
depends on the radius of the neutrinosphere, (£'), and is given by 
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Fig. 19. — The thick soHd Une and thick dashed hne represent the error in the matching of the terms (CfDlD\,0],) - 47rr(l + e + p/p)H 
and {C^D^D'^) respectively. All other matches are to machine precision. Therefore, the time integration of these energy rates must 
reproduce the drift in the total energy represented by a thick solid line in Fig. \\^\. The two monitored expressions do not allow a 
distinction of errors induced by approximations in the matching procedure in Eqs. I69i and 1841 from errors induced by the application 
of the adaptive grid to radiation quantities in the comoving frame (see discussion following Eg. 1461 *). In an attempt to disentangle these 
two contributions to the violation of energy conservation, we compare the cancellation in the terms {Cf D^D^Oj^) — 47rr(l + e + p/p)H 
and {C^D'f'D%) with a Lagran gian run with 400 zones and 2 angular bins (thin solid line and thin dashed line, respectively). The energy 
violation of the former term is greatly reduced on a fixed grid. 



and 



r J V RAE) 

The quantity /i is the cosine of the angle from the Umb of the neutrinosphere to the point at r, corrected for gravitational 
bending, and fiQ is that same angle as seen in the fluid frame. The net diffusivity in the diffusion equation is then set 
by the minimum of the two parts of the flux limiter, 

T) {E) = min (^i„trp {E) , ^gcom {E)) . 

In the following, we compare Bruenn's mgfld simulation with this flux limiter to the solution we have obtained by the 
complete solution of the transport equation. Beside of the different methods implemented for the radiation transport, 
the two codes also use different schemes to solve the hydrodynamics equations. While the code described in this paper 
adapts to spatial resolution requirements by continuously displacing zones, the MGFLD code adjusts the resolution 
by occasionally inserting or removing zones. It uses of order 150 zones at bounce and 250 zones half a second after 
bounce, while AGILE-boltztran works with 103 zones throughout. The energy resolution in mgfld is given by 20 
groups which adapt to the general relativistic redshift. In the Boltzmann solver, we used 12 fixed energy groups. 
We start with the comparison of a time slice at bounce in Figure (I20|) . The ordinates of graphs (a-h) display the 
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Fig. 20. — Comparison with MGFLD, 13 Mq model at bounce. 
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enclosed mass. At the end of core collapse, we find differences in the density profiles in graph (b) of 3% at the center, 
and less than 1% in the inner core. Differences of up to 80% around the shock front are due to a different resolution of 
the shock. Outside of the shock, we find a ~ 30% lower density in AGILE-BOLTZTRAN. The shock position, however, 
is in reasonable agreement as shown in the velocity profiles in graph (a). This is a consequence of the agreement in 
the thermodynamical state within the homologous inner core up to the sonic point at an enclosed mass of ~ 0.54 Mq. 
Differences in the infall velocities are of order 6%. In the diffusive domain, we find agreement in the entropies of order 
5% (graph (c)). Around an enclosed mass of 0.8 Mq, however, the deviations are systematic and of order 10%. At 
least part of it is due to a narrow choice of the spatial resolution in AGILE-BOLTZTRAN as discussed in section 14.51 
Immediately connected to the entropy profile is the temperature profile in graph (f). The temperature is very sensitive 
to entropy differences in the degenerate high density material. The electron fractions in graph (d) differ typically by 2%. 
The maximum deviation of 4% is found where the entropy deviation is largest. At lower densities than ~ 10^" g/cm^, 
the dominant source of neutrinos is provided by electron capture on heavy nuclei. As the mass fraction of heavy nuclei 
is not very temperature dependent, the entropy difference only affects the deleptonization by the dependence of the 
electron capture rates. AGILE-boltztran finds a slightly lower electron fraction in this regime. At densities exceeding 
^ 10^^ g/cm'^, however, standard input physics used to prematurely switch off electron captures on heavy nuclei in the 
oversimplified independent particle model iLanganke et al. 2003) such that electron captures on free protons became 
dominant in the simulations compared here. The lower entropy in the simulation with AGILE-boltztran leads to a 
smaller mass fraction of free protons and reduced deleptonization with respect to the MGFLD simulation. The slightly 
larger electron abundance between 0.65 and 0.8 Mq is a consequence. The agreement in the luminosities in graph 
(e) is of order 15%. The solution in the diffusive domain is smoother in AGILE-boltztran because the expanding 
zones of the adaptive grid equilibrate local discontinuities introduced by the equation of state. In contrast to these 
random fluctuations inside the shock front, we find a systematic deviation of 15% in the luminosities outside of the 
shock front. This difference might also be a consequence of the lower entropy in AGILE-BOLTZTRAN. The neutrino 
rms energies (graph (h)) in the inner core are determined by thermal equilibrium. There, the agreement is to 3%. 
Outside of the shock front, the deviation is initially more around 6%, before it becomes better again towards very 
large radii. The mean flux factor, the ratio of neutrino flux over neutrino density, F = H/ {cJ), in graph (g) shows 
differences of order 20% at this time. The smaller flux factor of the mgfld neutrinos at an enclosed mass of about 
0.8 Mq is probably due to increased isotropic neutrino emission at the higher entropy. It leads to slightly higher 
luminosities and higher rms energies, and is consistent with the smaller electron fraction. We should note here that it 
is very typical for such comparisons that each code presents an itself consistent picture of the dynamics such that it 
is sometimes almost impossible to isolate a single cause for a difference in the strongly coupled variables. The smaller 
flux factor in AGILE-boltztran outside 1.1 Mq stems from the limited angular resolution discussed in section l|4.2|l . 
It is systematic and visible in all the following figures. Its impact on the dynamics, however, is negligible because of 
the small coupling between the neutrino flux and matter at large radii. 

Figure H21|l shows a time slice just after the launch of the neutrino burst. This is the most dynamical phase for the 
radiation transport. The neutrino burst in graph (e) is evident in both codes as a propagating peak in the luminosity 
profile. The peak at ^ 2000 km radius is broader and by 25% smaller in AGILE-boltztran because of the numerical 
diffusion we have analyzed in section The mgfld burst changes its shape more slowly because the mgfld code 
reverts, in the transparent regime, to centered difference advection (rather than first order upwind) which is second 
order accurate in space (on a uniform spatial grid) . The luminosity in the region between the shock and the luminosity 
peak decays very rapidly. Differences of order 15% are due to a small time mismatch between the two solutions. 
Inside the shock front, the luminosity is still much noisier in the mgfld solution. This is also true for the entropy 
profile in graph (c), where local differences of 20% behind the shock or 50% at 400 km radius, disturb the otherwise 
nice agreement of order 5%. Note however, that the temperature profiles in graph (f) is less sensitive to this entropy 
variation. The temperature difference at 400 km radius is 12%. The agreement in the infall velocities in graph (a) 
has improved to generally 2% agreement with the exception of 9% around 400 km radius, i.e. the region with the 
entropy differences. The entropy deviations inside the shock front and at 400 km radius cause density differences of 
25% in these regions. The central density agrees to 1% and the agreement in the very distant layers is of order 3%. 
The electron fraction profile is affected by the numerical noise as well. Deviations, however, are only of order 2%, with 
local exceptions showing differences of 10%. the rms energies in graph (h) typically agree to 5%. The differences in 
the flux factors in graph (g) show systematic deviations of 10% outside of 100 km radius and random fluctuations of 
order 20% interior to it. The latter are probably a consequence of the variations in the density and entropy profiles. 

Most important for the success or failure of the supernova explosion in our spherically symmetric simulations is 
the time ~ 100 ms after bounce. Though, we have to state clearly that the general relativistic runs are not failing 
marginally. They are failing with officious insistence, and differences of the size we have described above are far from 
changing this, as we will show with the next two time slices. Figure H22I) shows the comparison of the two simulations 
at 100 ms after bounce, when the efficiency of the neutrino heating is close to maximum. The agreement of the data 
in this phase is remarkable. In this quasi-stationary phase, the position of the stalled shock is in accurate agreement 
as one can see in graph (a). We note however, that the important infall velocity in the heating region (e.g. at 100 km 
radius) is 10% larger in MGFLD than in AGILE-BOLTZTRAN. Since this infall velocity determines the time the infalling 
matter spends in the heating region before being accreted onto the protoneutron star, its value is as relevant as the 
heating rate. Better visible in the graph are the differences in the infall velocity which have only a size of 3%. Again, 
they reflect differences in the entropy profile in graph (c). Outside the shock the entropy deviations of order 20% 
seem large. However, they correspond to temperature variations of only 10% and are explicable by the assumption of 
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Fig. 21. — Comparison with MGFLD, 13 model, 10 ms after bounce. 




Fig. 22. — Comparison with MGFLD, 13 Mq model, 100 ms after bounce. 
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instantaneous silicon burning in AGILE-boltztran. mgfld uses a 9 species nuclear network to burn nuclei to iron 
until the temperature exceeds 5.106 x 10^ K, at which point the zone is flashed to NSE. Inside the shock front, the 
agreement in the entropies is good to 2% (except for some local fluctuations in the mgfld solution in the heating 
region). As the rise in entropy between the shock front and the location at 100 km radius is entirely due to neutrino 
heating, this is an encouraging result. In this time slice, we have replaced the temperature graph (f) with a graph 
of the more important heating rates. Plotted is the sum of the heating rates by electron neutrino absorption and 
electron antineutrino absorption. The rate is negative in the cooling region, representing cooling rates by the inverse 
processes, electron capture and positron capture. We find 7% difference in the heating region and 10% difference in 
the cooling region. The values in the mgfld simulation are larger. This is perfectly consistent with the larger electron 
antineutrino luminosity in graph (e) and the larger infall velocity behind the shock in graph (a). The electron fraction 
profiles in graph (d) agree to 2% accuracy with the exception of two local 6% deviations at 25 km and 45 km radius, 
where the electron fraction is as small as ^ 0.1. This close agreement in the entropy and electron fraction leads to 
a similar thermodynamical state of the fluid in the inner core. Good agreement in the nearly hydrostatic density 
profile in graph (b) is the consequence. Differences are 0.5% at the center and on average around 5% up to a radius 
of 400 km. Further out, differences stemming from the silicon/iron layer interface are visible. The good agreement in 
the density profiles facilitates the comparison of other density-dependent quantities. E.g. the neutrino rms energies 
in graph (h) agree to 2% in the innermost core, to 6% in the heating regions, and to 3% outside of the shock front 
in the neutrino signal. In the heating region, the rms energies are lower in AGILE-BOLTZTRAN. At larger radii, the 
electron neutrino energies are lower in MGFLD, while the fjL- and r-neutrino energies are lower in AGILE-boltztran. 
The electron antineutrino energies do not show a discernible tendency. The electron neutrino luminosities in graph 
(e) are 5% larger in AGILE-boltztran in the diffusive domain and 2% larger in mgfld in the heating region. The 
electron antineutrino luminosity is about 6% larger in AGILE-BOLTZTRAN in both domains. In the /i- and r-neutrino 
luminosities, deviations start at 20 km radius with 2%and linearly increase to 10% at the shock position. With respect 
to the neutrino signal, the electron flavor luminosities are larger in the simulation with MGFLD and the ^- and r- 
luminosities are larger in the simulation with AGILE-boltztran. This is just inverse to the relation between the rms 
energies. While, for example, an increase in the rms energy contributes only linearly towards a luminosity increase, 
the neutrino flux, however, is determined by the quadratically decreasing mean free path. Where the flux factors in 
graph (g) are below 0.5, they are typically 5% lower in mgfld. Between this domain and the shock front (including 
the heating region) the flux factors show on average no discernible deviation, but fluctuations of ±5% are evident. 
The flux factors in AGILE-BOLTZTRAN disagree with the asymptotic limit, F = 1, by 8% at distant radii due to the 
limitation in angular resolution as discussed in subsections 14 . 21 and 14 . 51 

Similar agreements and differences are found in later stages of the quasistationary phase, e.g. at 400 ms after 
bounce (Fig. (|23|l ). We find again hotter neutrinos between the neutrino sphere and the shock position in the MGFLD 
simulation and a larger diffusive flux at very high densities in AGILE-boltztran. On a longer time scale, the latter 
leads to a marginally larger deleptonization, visible in the innermost 7 km of the electron fraction profile. The electron 
fraction in AGILE-boltztran is now 2% smaller, whereas it was 2% larger early after bounce. A new difference of 
6% appears in the peak of the entropy profile, which might be due to differences in the hydrodynamics scheme for the 
solution of the shock jump conditions. We also point to the fact that the resolution of the shock in this stage is quite 
poor. The thin marks at the top of graph (a) represent the grid point locations in the AGILE-BOLTZTRAN simulation. 
Most of the grid points cluster around the forming density cliff between 20 km and 30 km radius. A high resolution in 
this domain is important for the radiation transport because most neutrinos are emitted from the cliff. However, this 
leaves less grid points to focus on the shock front. Graph (a) shows only a small increase of the grid point concentration 
at the shock front. The infall velocities inside the shock in graph (a) have increased to 10** km/s. At this speed, the 
infalling material is only very shortly exposed to neutrino heating. This is consistent with the observation that the 
entropy increase from the shock to the peak entropy is only marginal in both runs. Without unaccounted effects, a 
shock revival is impossible under these conditions. The most significant new differences appear in the luminosity and 
rms energy profiles, mgfld has now consistently smaller luminosities and higher rms energies. The differences are 
20% in the luminosities and unchanged 6% in the rms energies in the heating region. Also interesting to note is that 
the rms energies agree perfectly right outside the shock at 100 km radius, and develop deviations of 3% — 4% at far 
distances. There might be small differences in the implementation of the observer corrections. The flux factors in 
graph (g) are more difficult to compare, because the neutrinos no longer decouple at exactly the same radius. The 
flux factor in AGILE-BOLTZTRAN is still numerically limited in its asymptote, while the flux factor in mgfld exhibits 
somewhat peculiar fluctuations in the heating region. 

Overall, in our general relativistic simulations of stellar core collapse and postbounce evolution in spherically symmet- 
ric space-time, we find quite exactly the same physical evolution in both runs. One implementing the full Boltzmann 
transport equation, and the other based on the multi-group flux-limited diffusion approximation with a new sophisti- 
cated flux limiter. This is especially true in all phases of stationary radiation transport. There are modest quantitative 
differences in these simulations of failed supernovae. However, except for some obvious cases, like e.g. the limited 
asymptotic flux factor in AGILE-BOLTZTRAN, or the numerical fluctuations that the equation of state induces to purely 
Lagrangian hydrodynamics in MGFLD, these quantitative differences cannot uniquely be traced back to principal weak- 
nesses or strengths of cither method. They may rather be due to method-specific details in the implementation of two 
very different methods to solve a complex time-dependent physical problem. The advantage of the Boltzmann solution 
is that it is complete, and, beside of numerical resolution, does not need to justify or verify any basic assumptions. 
The advantage of mgfld is the transparency of the approach, and that it produces accurate results with much higher 
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efficiency in the application we have tested here. 

5. CONCLUSION 

AGILE-BOLTZTRAN directly solves the general relativistic Boltzmann transport equation for the specific neutrino 
distribution function in spherically symmetric space-time. In combination with general relativistic hydrodynamics, 
we present a consistent finite differencing for fermion radiation hydrodynamics in spherical symmetry. We test our 
code in the context of stellar core collapse and postbounce evolution of a smaller 13 Mq, and a more massive 40 Mq 
progenitor star. They embrace the progenitor mass range expected to produce explosions after core collapse. But no 
explosions are obtained in our simulations. We analyze the regions of neutrino emission by a statistical description 
allowing the detailed presentation of the specific production reactions and energies of escaping neutrinos. 

The Boltzmann transport equation plays a fundamental role in kinetic gas theory. One may derive from it the 
diffusion equation, or — if one drives the limit further to vanishing transport — the hydrodynamics equations. Thermo- 
dynamical quantities are given by the expectation values of various operators on the fundamental particle distribution 
functions. It is quite unique to numerically solve a single equation for the fundamental particle distribution functions 
in an object of astronomical size! However, this basic approach becomes intricate in the detailed finite difference 
representation. In theory, the solution of the Boltzmann equation guarantees the correct evolution of expectation 
values in the energy conservation equation or the diffusion limit. Also the equation of state for radiation, p oc p/3, will 
emerge correctly. This is not automatically the case in a numerically solved Boltzmann equation. Especially in the 
astrophysical application, a distant observer will rather observe macroscopic phenomena than the microscopic local 
state of the matter and the radiation field. Therefore we compose the finite difference representation of the micro- 
scopic physics such that important macroscopic quantities are accurately met. This optimizes the accuracy of the 
simulations at moderate resolution settings. However, it would be more straightforward to comply with this ambition 
if there would not be an additional complication. Physically, the Boltzmann transport equation is a trivial equation, 
describing free particle propagation on geodesies between independent collisions (even more so in spherical symmetry, 
where only two independent spatial degrees of freedom have to be considered). With respect to the collisions, the 
supernova has been investigated as a site of interesting nuclear physics and neutrino-matter interactions. The cross 
sections depend on the local neutrino energies and scattering angles. They are most conveniently evaluated in the 
fluid rest frame and most easily described in combination with comoving coordinates. The latter are also favored by 
the initial core collapse where the relevant computational domain collapses with the matter from several thousands to 
some hundreds of kilometer radius, before it hopefully expands again with the explosion. The comoving coordinates 
convert the simple statement of the Boltzmann equation, that the derivative of the invariant particle distribution 
function along the phase flow vanishes between collisions, into a complicated sum of partial derivatives along the 
comoving coordinates. The finite differencing of all comoving frame correction terms have to be adjusted in a mutually 
dependent way to correctly transform to conservation equations for the observer in a laboratory frame. We construct 
number conservation to be guaranteed. With respect to energy conservation, we match the largest 0{v/c) terms to 
machine precision, while higher order gravitational terms are matched approximately, guided by specific limiting cases. 
The adaptive grid applied to the neutrino distribution function in the comoving frame introduces additional deviations 
from perfect energy conservation in the laboratory frame. In the worst case, i.e. in the evolution of the very massive 
40 M0 progenitor star, all these deviations together reach the order of a supernova explosion energy 10^^ erg) 
at the time the protoneutron star collapses to a black hole. However, when the conditions are most favorable for an 
explosion (we should rather say least unfavorable), the numerical energy gain has not yet exceeded 10^° erg and we 
are confident that the remaining energy drift does not influence our physical conclusions. We have also tested the 
implementation of the gravitational redshift and bending, the angular advection in the free streaming regime, and 
the propagation of the neutrino burst. Quite generally, we find that the choice of first order upwind differencing is 
the limiting factor for the accuracy in our numerical solution. In space, it causes artificial pulse spreading. In the 
advection terms of the adaptive grid, it causes artificial diffusion. In momentum space, it causes angular diffusion in 
the evolution of the particle distribution function. Although these limitations can be reduced by the choice of higher 
resolution, more elaborate advection schemes could generally be beneficial for the accuracy of the simulations. With an 
investigation of the resolution dependence of the results in each phase space dimension, however, we demonstrate that 
our physical results are sufficiently converged in the supernova application. The undesired effects of low order upwind 
differencing predominantly affect the regions far from the ncutrinospheres, where the radiation field is decoupled from 
the matter such that an adverse influence on the dynamics of the model can be excluded. We derived the moments 
entering the energy conservation equation. Moreover, we algebraically derived the flnite difference representation of 
the diffusive limit, the nonrelativistic limit, and the radial dependence of the radiation quantities in a stationary state 
free streaming radiation field. These algebraic gymnastics led to a better understanding of the internal mechanisms 
in our finite difference representation and, by comparison with the code results, enhanced the confldence that the 
implementation exactly corresponds to the programmer's intention. In order to perform simulations over a second 
or more, the radiation transport and hydrodynamics have been implemented with implicit flnite differencing. This 
allows reasonably sized time steps during the neutrino heating phase. Finally we compared the evolution of a 13 Mq 
star with an independently implemented general relativistic supernova code which applies the multi-group flux-limited 
diffusion approximation with a recently developed ffux limiter (Brucnn 2002). We find significant agreement that is 
especially impressive in the neutrino heating phase. Other phases show regions with moderate deviations, but those 
are far from affecting any physical conclusions. 

With the exponentially growing power of computer hardware, computational astrophysics may reach a comparable 
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status in astronomical, nuclear and particle physics research as terrestrial experiments have received in the previous 
century. We may sort scientific progress into three phases: First, the astronomical observation or a new theoretical 
concept lead to a primary idea about the basic physics involved in an event. At this stage, order of magnitude 
estimates are made and compared to the qualitative observational data. The second step involves plausibility studies 
of the suggested scenario, which may include detailed numerical simulations and comparisons to observations. Often, 
however, approximations have to be used that are dictated by technical limitations in the newly developing field. Only 
in the simplest cases it can be shown that all assumptions are based on undeniably justifiable physical considerations. 
This exciting phase of discovery and controversy should be followed by more rigorous numerical simulations in a 
third step, when technical limitations fade and remaining approximations become quantifiable. This phase has to 
be accompanied by code documentation that allows independent researchers to analyze and eventually reproduce 
the numerical results, similar to the way the value of experimental data is enhanced by a detailed description of the 
reproduceable experiment. The recent emergence of functional Boltzmann solvers leads the radiation hydrodynamics in 
spherically symmetric supernova models into this third phase. The third phase, by principle, should be accompanied by 
solid agreement in the numerical solution found with independent implementations of the same physical ingredients. 
In order to perform such comparisons, we appreciate the situation that a basic set of standard nuclear and weak 
interaction phy sics has been defined and used for over 15 years (e.g. ijTubbs fc Schrarmn. 1975: Schinder fc Shapirol 
119821: iBruenn i lil985'l'l. In a comparison with the second phase — which used the former standard scheme for neutrino 
transport in supernovae — multi-group flux-limited diffusion, we hope to increase the confidence in our results. On 
the other hand, we support the validity of the MGFLD approximation in spherical symmetry if the flux limiter is 
chosen carefully. As satisfying as these close results are, they should not be mistaken as an example of convergence 
to agreement in the sense of the third phase. The third phase is only reached if similar agreement is found between 
independent groups and other documented codes that solve physically complete neut rino transport equations with 
general relativistic effects, e.g. the variable Eddington factor method implemented bv iRampp & Janka I l)2fl02D . We 
hope that spherically symmetric simulations with neutrino radiation hydrodynamics have reached phase three and 
that accurate and detailed neutrino information continues to be useful for the exploration and improvement of the 
local microscopic input physics. 
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APPENDIX 

0(y/C) LIMIT OF THE FINITE DIFFERENCE REPRESENTATION 

The 0{v/c) Boltzmann equation has been derived by Castor ijCastor Il972() . One can also obtain it by the elimination 
of higher order terms from Eq. 1)15(1 . It is enough to set a = T = const. = 1 and to replace u by the nonrelativistic 
velocity v. The conser vation properties of the O jvlc) Boltzmann equation become apparent when we take its energy 
and angular moments ((Mihalas fc Mihalas II198J) : 

9 J d r 2 ttI ^ / t f dlnp 2u\ ^ 

-E^dEdfi + / xFE^dEdfi = 0, 
P J 



9H d . 1 , ^ fd\np 2v 



r 

+ j xFE^dEpdp = 0. (Al) 

As we construct the specific radiation energy in the laboratory frame, J + vH , we keep all terms that arise from the 
0(v/c) Boltzmann equation and obtain, in analogy to Eq. ((27|l . an almost conservative evolution equation for the 
radiation moments, 

= — ( J + wff ) + — Unr^p {vK + H)] 

■^E^dEd^i + J xFE^dEdjjL + v j xFE^dEndn 

1 1 

(47rr^pv)H. (A2) 



Anr'^p dt 
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Although terms of higher order than (w/c), especially the gravitational terms, require consideration in realistic su- 
pernova simulations, the 0{v/c) limit may be useful for various explorative studies. We give below a finite difference 
representation of the 0(w/c) Boltzmann equation by simply substituting a = F = 1 in the corresponding general 
relativistic expressions: 

F, I, -F; ■/ 1 
t' .-j'-k' ^ r .f ,k fin , 1 r 1 1-1 

^* ^ — — dt + d^, ["^+1^+1./-^' - (A3) 
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Pi' 

The expressions basically reduce to the finite difference representation described in (jMezzacanna fc Bruenn lll993a() . 
except for the adaptive grid extension, the new code flow outlined in section 13.51 the improved choice of transport 
coefhcients (3i^k' and 7i',fc' in Eq. (|56|l . the new discretization of angular aberration, and the matched finite difference 
representation of Ai' ^k' and Bi'_ji^k' , 

47r p^' 

Ai'^k' = " (r-f+i {vi+2 - Wi+i) /3i+i,fc' + {vi+i - Vi) (1 - /3i,fc')) (A9) 

B^'.3',k' = _ 1 '^^ [7»',fc'Cj+i {Pj'+i - H') + (1 - 7i',fc') Ci {^^f - Mj'-i)] ■ (AlO) 
The angular difference coefficients, Cj, are defined in Eq. (|63|l . 

ATTENUATION FACTORS FOR THE PRESENTATION OF INTERACTION RATES 

Whenever one composes a graph for the discussion of weak interaction rates in the supernova environment, one has 
to circumvent the inherently large scale differences of the rates at different locations in the star. In a logarithmic 
presentation, many details are hidden and differences between absorption and emission are difficult to appreciate. One 
can focus on the neutrinospheres and investigate the rates of interest under the corresponding conditions. However, 
the definition of the neutrinosphere is only based on the opacity and does not account for large emissivities outside the 
neutrino sphere. One may miss important sources of the total neutrino luminosity and fail in the explanation of the 
spectra. Moreover, it is common to average the extremely energy-dependent location where a given optical depth is 
reached to one single neutrinosphere — an even more problematic concept. In this appendix, we motivate a convenient 
presentation of interaction rates according to their relevance to the total luminosities. The approach aims to produce 
intuitively accessible figures with information about where the neutrinos come from, which reactions contribute to the 
total luminosity, how they locally compare to other reactions, and how the neutrino spectra are formed. 

We derive auxiliary attenuation factors in terms of a staggered grid with zone edge indices i and zone center indices 
i' = i -\- 1/2. We start with a conserved luminosity Li with respect to spheres around the symmetry center. We may 
choose the neutrino number luminosity in units of particles per second or the neutrino energy luminosity in ergs per 
second. Each zone may have a number or energy source, erriii , and a number or energy sink, ahi' , in number or ergs 
per gram and second. The luminosity is then recursively defined by its central value Lq = and 

Li+i = L, + {errii' - ah') dai', (Bl) 

where daii denotes the rest mass contained in the mass shell i'. The entering luminosity, Li, and the source, emi'dai/, 
are subject to absorption in this shell. We define an attenuation factor, Xi' < 1, which accounts for the reduction of 
these quantities, 

Li+i = Xi> (Li + emi>dai,) . (B2) 
From Eq. ljBl(l and ljB2|l , we can readily isolate Xi' : 

Li+i Li^i 

Xi' —— . 

Li + emi'dai' Li^i -\- abi' dae 
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In the rare cases of a negative Xii , we set Xi' to zero. If we calculate a mean free path, A^', from all reactions that may 
act as a sink, the absorption rate can also be derived from the local specific number or energy density, Ji' , according 
to abi' = cJi'/Ai'. The attenuation factor Xi' can then be expressed by physical quantities that are well accessible in 
a numerical evolution of radiation hydrodynamics, 

(B3) 



On the other hand, we derive from Eq. ljB2p by two recursive self-substitutions, 

Li+i = Xi' {xi'-i {xi'^2 {Li-2 + eTOi/_2dai'_2) + eruif ^idai> ^i) + emi^dai') 
The continuation to = and a rearrangement of the terms leads to 



n / n 



i—1 / 

This equation describes how the total number or energy luminosity at a radius r„+i is composed by contributions of 
distributed emissivities in the star. The attenuation coefficients 

n n J 

l=^ i=t ^'+1 + -jfrdai' 

suppress irrelevant sources that are subject to large rea bsorp tion. It is instructive to go a step further. We introduce 
the flux factor hi' = Li+i/ (^Anr^cpJi'^ and rewrite Eq. IjBSp as 

dr,: ^ ~^ 



in order to extract the logarithmic derivative of ^, 

^n+l,i' + l — £,n+lA' dri' 



This is a finite difference representation of the equation 

rfln^ _ J_ 
dr hX 

with the solution ^ = exp ^— J (/lA) ^dr^. It is the familiar attenuation exp(— r), if the radius-dependent fiux 
factors are properly accounted for in the evaluation of the optical depth, t = J (hX) ^ dr. However, we evaluate the 



attenuation coefficients according to Eq. I|B5|) where the finite differencing is consistent with the conse rvation laws 
in our implementation of the Boltzmann equation. As a consistency check, we compare in Fig . (jB24|l the original 
luminosities Ln+i from the simulation with the luminosities reconstructed according to Eq. (|B4|) . 
We can use the attenuation coefficients in many convenient ways. A graph showing 

dtti' 

ffri+i [ri') = t,n+i.i'emi'- — 
art' 

as a function of radius, r^/, visualizes the contribution of each region in the star towards the total luminosity (represented 
by the area under the graph) at radius r„+i. In a comparison of different neutrino species, the graphs illustrate the 
decoupling at different radii. If the gravitational well is not too deep, the neutrino energy is approximately a constant 
of motion and the neutrino in different energy groups can be treated like different species. Moreover, instead of 
considering total emissivities and opacities, one can disentangle them into different reactions which we will enumerate 
with a superscript £. A figure with graphs showing 

gi+i{n') = ^n+iyemi^ (B6) 
dr,> 

as a function of radius r^' would visualize the contribution of reaction £ to the total luminosity at radius rn+i by 
the enclosed area under the line gfi^i. The graph is automatically scaled such that the most important reactions 
for the total luminosity are presented most prominently. If a reaction conse rves the analyzed quantity we ar e fre e 
to include or not include it in the evaluation of the mean free path in Eq. (|B5|I and the emissivities in Eq. ljB4|l . 
Scattering, for example, does not change the neutrino number. If we include scattering as a reaction in the analysis 
of the number luminosity, the attenuation coefficients would indicate the probability to escape from a given location 
without any further scattering. If we do not include scattering, the attenuation coefficients indicate the (in many cases 
larger) probability to escape from a given location without being absorbed on the way out. The number luminosity 
is the same in both cases because scattering conserves the number of propagating neutrinos. If we work with the 
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Fig. B24. — The thin lines show the luminosity profiles in the evolution of the 13 Mq progenitor star at 100 ms after bounce. The 
solid line represents the electron neutrino luminosity, the dashed line the electron antineutrino lumino sity, and the dash-dotted line the fi- 
and T- neutrino luminosity. The thick lines show the luminosity profiles evaluated according to Eq. 1B41 . i.e. based on attenuated local 
emissivities. The reconstruction obviously misses the observer corrections in regions with low interaction rates. Otherwise, it reproduces 
the original luminosities sufficiently well to be accurate in the analysis of the formation of the neutrino spectra. 



energy luminosity instead of the number luminosity, we have to include neutrino-electron scattering because this 
reaction affects the energy of escaping neutrinos after their production. We simply decompose a neutrino-electron 
scattering reaction into a neutrino absorption at the incoming neutrino energy and a neutrino pro duc tion a t th e 
outgoing neutrino energy. These terms are then included in the opacities and emissivities in Eqs. (|B5p and IIB4|l . 
In this case, the attenuation coefficients indicate the probability to escape from a given location without an energy- 
changing reaction. We found this choice to be the most interesting for the analysis of the formation of the neutrino 
spectra. The omitted isoenergetic scattering reactions are reflected in the transport spheres which are easily displayed 
as a complement. Finally, we mention the potential use of the attenuation coefficients for statistical evaluations, e.g. 
for the average radius of neutrino emission, 

(r) = ^JTn' , (B7) 

In contrast to the classical definition of the neutrinosphere, the quantity (r) accounts for regions with high emissivities 
in transparent regimes. Many other statistical informations at the origin of the neutrino emission may be obtained 
analogously. 
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